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I.  INTRODUCTION 


A.  THE  CONTROL  SYSTEM  DESIGN  PROCESS 

Since  the  beginning  of  time,  man  has  sought  ways  to  control  the  laws  of  nature. 
From  the  simple  float  regulator  developed  by  the  Greeks  in  300  B.C.  [Ref.  1:  p  3],  to 
the  amazingly  complex  space  shuttle  of  the  1980's,  control  systems  span  the  range  of 
mankind's  efforts  to  govern  his  surroundings.  The  challenge  for  a  control  systems 
engineer  is  to  use  his  knowledge,  skill,  judgment,  and  experience  to  systematically 
develop  a  solution  to  any  of  a  number  of  different  types  of  control  problems.  There  is 
seldom  only  one  right  answer  to  a  control  problem.  In  general,  there  may  be  several 
alternate  solutions  to  the  same  problem  and  the  final  product  will  probably  be  a 
compromise  between  them.  It  remains  the  responsibility  of  the  engineer  to  choose  the 
"best"  solution  that  meets  the  performance  criteria  specified  by  the  user.  So,  how  does 
the  engineer  know  where  to  begin  when  he  is  given  a  set  of  performance  criteria  for  a 
system  ?  There  is  no  set  procedure  carved  in  stone.  There  are,  however,  a  few  broad 
guidelines  that  give  the  engineer  a  rough  idea  of  the  tasks  which  need  to  be 
accomplished  in  his  quest  to  design  an  effective  control  system.  These  milestones  are 

as  follows : 

% 

1.  Define  the  system. 

2.  'Specify  the  desired  performance  of  the  system. 

3.  Identify  the  constraints  under  which  the  system  must  operate. 

4.  Translate  the  information  from  milestones  1,  2,  and  3  into  a  mathematical 
model  that  can  be  simulated  on  the  computer. 

5.  Use  the  available  tools  to  develop  a  control  system  which  satisfies  the 
performance  specifications. 

6.  Evaluate  the  control  system  design  using  computer  simulations. 

7.  Modify  the  design  as  required  to  better  suit  the  application. 

8.  Incorporate  the  control  design  in  a  prototype  system  to  test  the  ability  of  the 
system  to  tolerate  real  world  non*linearities  and  non-ideal  conditions. 

9.  Modify  and  optimize  the  design  until  a  satisfactory  control  system  is  realized. 
The  first  milestone  listed  above  is  not  always  as  simple  as  it  appears  to  be.  In 

fact,  defining  the  limits  of  the  system  may  possibly  be  the  most  difllcult  phase  of  the 
design.  If  the  engineer  docs  not  expend  considerable  effort  in  the  definition  of  the 
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problem  which  he  is  to  solve,  he  might  end  up  solving  the  wrong  problem.  The 
engineer  needs  to  include  enough  parameters  in  his  system  to  accurately  model  the  true 
system  without  becoming  overburdened  computationally.  This  is  as  much  an  art  as  it 
is  a  science.  The  engineer  will  probably  need  to  make  simplifying  assumptions  and 
approximations  which  tend  to  widen  the  gap  between  the  performance  of  the  model 
and  the  performance  of  the  real  world  system. 

Defining  the  desired  performance  of  the  system  may  or  may  not  be  left  to  the 
discretion  of  the  engineer.  If  a  strict  set  of  specifications  is  handed  to  him,  then  the 
engineer  has  little  choice  but  to  satisfy  those  specifications  or  bei  able  to  defend  his 
claim  that  they  can  not  be  satisfied.  On  the  other  hand,  there  may  be  considerable 
leeway  for  the  engineer  to  make  sweeping  changes  in  the  control  system  and  still  satisfy 
the  required  specifications.  Several  tools  are  available  to  measure  the  performance  of  a 
control  system.  The  classical  design  engineer  holds  fast  to  such  measures  as  the  gain 
margin,  phase  margin,  root  locations  in  the  S  plane  or  Z  plane,  and  bandwidth.  The 
advent  of  optimal  control  techniques  has  placed  emphasis  on  the  minimization  of  some 
cost  function  as  a  means  of  measuring  system  performance.  All  of  these  techniques 
have  their  place  in  the  realm  of  control  system  design  and  it  is  the  mark  of  a  successful 
designer  that  he  can  incorporate  any  or  all  of  the  tools  when  the  situation  dictates. 

The  third  milestone  is  an  important  yet  often  overlooked  element  of  the  design 
process.  Constraints  on  the  system  may  include  any  or  all  of  the  following  : 

1.  Monetary  cost 

2.  Admissible  control  inputs 

a.  Saturation  limits 

b.  Observability  of  the  parameters  required  for  control 

3.  Physical  limitations 

a.  Size 

b.  Weight 

c.  Minimum  and  /  or  maximum  velocities,  accelerations,  etc. 

d.  Imtial  conditions 

e.  Final  conditions 

f.  Sampling  rate  and  processor  speed  for  digitally  controlled  systems 

The  mathematical  model  is  the  link  between  the  real  world  system  and  the  design 
tools  which  the  engineer  has  at  his  disposal.  In  general,  most  physical  systems  which 
need  to  be  modelled  can  be  represented  by  some  set  of  differential  equations.  The 
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physical  laws  of  nature  lend  themselves  nicely  to  approximation  by  linear  ordinary 
difTerential  equations  with  constant  coefficients.  Non-linearities  and  random  processes 
are  also  quite  prevalent  in  many  systems  and  the  effects  of  these  phenomenon  can 
greatly  complicate  the  engineer's  effort  to  model  a  system.  That  is  why  he  must  have 
the  expertise  and  experience  to  know  how  to  make  assumptions  and  approximations 
which  simplify  the  problem  at  hand  to  a  point  where  he  can  use  the  available  design 
tools. 

The  next  step  is  to  use  design  tools  to  develop  a  control  system  which  satisfies 
the  desired  performance  specifications.  It  is  at  this  point  that  two  schools  of  tlTought 
begin  to  emerge.  The  classical  school  of  thought  focuses  on  such  design  tools  as  the 
Root  Locus  Plot,  Nyquist  Plot,  Bode  Diagram,  and  Function  Minimization.  The  more 
daring  school  of  thought  centers  its  attention  on  the  maximum  principle  of  Pon  tryagin 
and  the  method  of  dynamic  programming  developed  by  Bellman.  The  advent  of  digital 
computers  in  the  mid  1950's  made  these  more  powerful  design  tools  realizable  since  the 
amount  of  computation  required  by  them  was  prohibitive  if  not  impossible  to  do  by 
hand,  In  any  event,  the  design  engineer  has  a  myriad,  of  tools  from  which  to  choose. 

Once  the  design  of  the  controller  is  accomplished,  the  next  step  Is  to  integrate  the 
control  system  with  the  system  model  and  then  simulate  the  entire  system  to  evaluate 
its  performance.  A  very  useful  method  to  perform  this  evaluation  is  to  study  the  time 
response  of  the  output  variables  of  the  system.  Such  parameters  as  rise  time,  peak 
overshoot,  and  settling  time  are  typical  values  to  be  noted.  The  digital  computer  once 
again  is  tr  very  useful  means  of  obtaining  such  information  rapidly. 

Even  the  best  of  control  designers  is  not  apt  to  hit  a  buUseye  on  his  first  shot. 
Control  system  design  theory  does  not  guarantee  success  on  every  try.  The  method  of 
trial  and  error  is  one  with  which  all  control  engineers  are  familiar.  Modifications  to 
the  control  system  are  inevitable. 

Once  the  controller  design  is  proven  in  simulation  studies,  it  is  time  to  test  it  out 
on  a  prototype  or  small  scale  model  of  the  actual  system.  This  phase  of  design  can 
become  costly  if  the  designer  has  not  thoroughly  tested  his  controller  on  the  computer 
first.  !n  this  phase  of  design,  the  non-linearities  and  random  effects  which  were 
ignored  or  approximated  during  the  modelling  phase  become  significant  factors  once 
again.  Conditions  which  were  assumed  to  be  ideal  in  the  model  now  become  non-ideal. 
The  evaluation  process  begins  all  over  again  as  these  new  disturbances  change  the 
performance  of  the  system. 


13 


The  next  step  is  obvious.  After  evaluation  of  the  controller  in  a  real  world 
system,  the  need  for  further  modifications  is  again  probable  if  not  inevitable.  Changes 
must  be  made  until  a  satisfactory  controller  has  been  designed  which  meets  the  desired 
specifications. 

B.  AROD 

It  is  the  goal  of  this  thesis  to  complete  the  first  seven  steps  of  the  nine  step 

design  process  discussed  in  the  previous  section.  The  system  chosen  for  this  endeavor 

is  an  airborne  remotely  piloted  vehicle  (RPV)  called  AROD.  The  acronym  stands  for 
<  .  • 

Airborne  Remotely  Operated  Device.  The  United  States  Marine  Corps  initiated  work 
on  AROD  early  in  1986  and  is  attempting  to  introduce  the  vehicle  into  the  operational 
Fleet  Marine  Force  during  fiscal  1987.  AROD  is  a  slow,  low-flying  ducted  fan  vehicle 
pow'cred  by  a  vertically  mounted,  two  cycle,  two  cylinder  gasoline  engine  which  drives 
a  three- bladed  propeller.  See  Figure  1.1. 

The  vehicle  is  38  inches  tall,  32  inches  in  diameter,  and  has  a  nominal  weight  of 
85  pounds.  It  presently  has  a  payload  capacity  limited  to  a  miniature  television  camera 
and  a  canister  of  fiber-optic  cable.  The  Marine  Corps  plans  to  use  AROD  for  short 
range  reconnaissance  and  over-the-hill  spy  in  the  sky  surveillance.  The  fiber-optic 
cable  provides  two-way  communication  between  the  vehicle  and  the  ground  based 
operator.  The  uplink  communication  will  consist  of  control  commands  while  the 
downlink  will  provide  real  time  surveillance  and  on-board  status  information. 

The  primary  flight  mode  is  low  altitude  hovering  with  the  axis  of  the  spinning 
propeller  oriented  perpendicular  to  the  surface  of  the  earth.  In  order  to  translate 
horizontally  across  the  earth's  surface,  the  entire  vehicle  must  be  tilted  towards  the 
direction  of  intended  mevenrient.  The  mechanism  by  which  AROD  is  tilted  for  such 
translation  consists  of  four  control  vanes  located  in  the  airflow  downstream  of  the 
propeller  wash.  One  pair  of  control  vanes  is  designated  as  the  rudder.  The  other  pair 
is  designated  as  the  elevator  and  is  oriented  such  that  its  axis  of  rotation  is 
perpendicular  to  the  axis  of  rotation  of  the  rudder  pair.  See  Figure  1.2.  All  four  fins 
assume  dual  responsibility  for  aerodynamic  control  in  that  they  also  serve  as  ailerons 
for  AROD.  The  control  vanes  are  actuated  by  model  airplane  servos.  These  servos 
are  limited  to  a  maximum  deflection  of  ±  30®  and  a  maximum  angle  rate  of  50«/sec. 
A  maximum  translational  velocity  of  30  knots  (  34  mph  )  in  a  no  wind  condition  is 
desired.  The  translational  velocity  of  AROD  is  proportional  to  the  tilt  angle  created 
by  the  rudder  and  elevator  control  vanes. 


Figure  1.1  Schematic  Drawing  of  AROD 


Figure  1.2  AROD  Control  Vanes 


Vertical  flight  control  of  AROD  is  accomplished  through  a  throttle  controller 
which  incorporates  the  same  type  of  model  airplane  servo  used  to  actuate  the 
aerodynamic  control  vanes.  The  throttle  controller  increases  or  decreases  the  engine 
RPM  as  required  to  raise  or  lower  AROD  vertically. 

C.  THE  PROBLEM 

AROD  is  interesting  from  the  standpoint  of  a  control  system  design  for  several 
reasons.  Most  significant  is  the  phenomenon  of  cross>coupled  dynamics  between  the 
pitch  and  yaw  subsystems.  In  addition,  AROD  is  a  Multi-Input  Multi-Output  system. 
These  topics  arc  briefly  discussed  in  the  following  sections. 


1.  Gyroscopic  Coupling 

AROD's  propeller  has  a  spin  velocity  of  7200  RPM  in  the  hovering  condition. 
This  creates  a  large  angular  momentum  vector  along  the  spin  axis  of  the  propeller. 
Thus,  AROD  can  be  thought  of  as  a  large  gyroscope  with  its  angular  momentum 
vector  oriented  perpendicular  to  the  surface  of  the  earth.  Consider  a  cylindrical  rotor 


Figure  1.3  Spinning  Rotor  Orientation 

spinning  about  the  x-axis  with  an  angular  spin  velocity  (O.  See  Figure  1.3.  Let  the 
rotor  be  of  mass,  m,  with  moments  of  inertia  h  respective 

axes.  These  moments  are  defined  in  the  three  equations  below. 


JJJ  S(yJ  +  zh  dV 

(1.1) 

f dV 

(1.2) 

JJJ  8(x^  +  y*)  dV 

(1.3) 
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In  these  equations,  d  is  the  density  of  the  rotor  and  dV  is  an  incremental  volume 
element  of  the  rotor.  In  the  case  of  AROD,  it  is  assumed  that  the  propeller  is  spinning 
with  sufficient  angular  velocity  that  its  dynamics  can  be  approximated  by  the  dynamics 
of  a  cylindrical  disk  having  the  same  mass,  m,  radius,  r,  and  thickness,  h.  The 
moments  of  inertia  of  the  propeller  then  reduce  to  the  following  : 


mr2 

2 


(1-4) 


ly- 


mOr^  +  h^) 
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(1.5) 


mOr^  +  h^) 
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(1.6) 


Notice  that  the  moments  of  inertia  about  the  y-axis  and  z-axis  are  equal  to  each  other 
due  to  the  symmetry  of  the  problem. 

The  angular  velocity,  <d  ,  of  the  rotor  induces  an  angular  momentum  vector, 
H^,  defined  by  the  following  equation. 


W*  •  >x  <“ 


(1-7) 


If  a  coupled  pair  of  forces,  F,  directed  parallel  to  the  z-axis  is  applied  to  the  the  spin 
axis  of  the  rotor  at  a  distance,  d,  from  the  rotor's  center  of  gravity,  as  in  Figure  1.4, 
then  a  torque,  M  ,  results-  The  torque  is  given  by 


M  -  F  X  d 


(1.8) 


If  the  rotor  were  not  spinning  wi^h  angular  velocity,  co  ,  then  this  applied  torque  would 
result  in  rotation  of  the  rotor  about  the  y-axis.  The  angular  momentum  of  the 
spinning  rotor,  however,  results  in  quite  a  different  response  to  the  applied  torque. 
According  to  Newton's  Second  Law,  an  external  force  applied  to  the  center  of  gravity 
of  a  rigid  body  results  in  a  change  in  the  velocity  of  that  body.  A  corresponding 
change  in  the  body's  momentum  also  results.  The  changes  in  velocity  and  momentum 
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—  Figure  1.4  Spinning  Rotor  with  a  Force  Couple  Applied 

are  in  the  direction  of  the  applied  force.  This  law  extends  into  the  realm  of  angular 
forces,  or  moments,  and  angular  velocities.  In  short,  an  applied  moment,  M  ,  results 
in  a  change  in  angular  momentum. 


M 


dl^  (Q 
dt 


(1.9) 


Notice  that  the  change  in  angular  momentum  is  in  the  same  direction  as  the  applied 
moment.  This  is  the  key  to  gyroscopic  precession.  By  vector  addition  of  and  M  , 
it  can  be  seen  that  a  new  angular  momentum  vextor,  ,  results.  The  new  angular 
momentum  is  displaced  by  an  angle,  y,  from  the  initial  angular  momentum  vector,  H^. 
This  angular  movement  is  called  precession  and  it  occurs  at  a  precession  rate,  r, 
oriented  as  shown  in  Figure  1.5  and  defined  by  Equation  1.10. 
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Figure  1.5  Resultant  Angular  Momentum  With  an  Applied  Torque 


(1.10) 


It  can  be  shown  that  H^,  M  ,  and  r  are  always  mutually  perpendicular  to  each  other 
[Ref.  2:  p.  335],  and  that  these  vectors  are  related  by  the  expression 

M-I,  r  (l.ll) 

The  handy  mnemonic  for  this  relationship  is  that  "spin  follows  torque".  In  this 
example,  the  spin  vector,  ,  that  results  from  the  applied  torque,  M  ,,is  rotated  in 
the  direction  of  the  torque.  Notice  that  the  magnitude  of  the  precession  velocity  is 
directly  proportional  to  the  magnitude  of  the  applied  torque. 


In  the  case  of  designing  a  control  system  for  AROO,  the  preceding 
development  is  quite  important.  Recall  that  the  prop'jller  spins  ct  a  noidnal  velocity 
of  7200  RPM.  The  angular  momentum  of  this  high  speed  rotor  bas  significant  effect 
on  the  fight  dynamics  of  AROD.  In  order  to  change  the  orientation  of  this  angular 
momentum,  as  is  required  to  accomplish  translational  flight,  a  considerable  torque 
must  be  applied.  This  torque  is  produced  by  the  four  control  vanes  iocated  in  the 
propeller  downwash.  Note  that  the  gyroscopic  nature  of  AROD  introduces  cross- 
coupling  of  the  pitch  and  yaw  dynamics.  For  instance,  when  a  pitching  torque  is 
commanded  about  the  y-axis  via  the  elevator  vanes,  AROD  must  first  undergo  an 
initial  yawing  motion  about  the  z-axis.  Similarly,  a  yaw  command  from  the  rudder 
vanes  results  in  an  initial  pitching  motion  about  the  y-axis.  These  cross-coupled 
dynamics  must  be  considered  by  the  engineer  when  designing  the  controller  for  the 
elevator  and  rudder  vanes. 

2.  Multiple  Control  Loops 

Another  feature  which  makes  AROD  interesting  for  the  control  engineer  is  the 
Multi-Input  Multi-Output  (MIMO)  nature  of  its  dynamics.  There  are  basically  four 
subsystems  which  need  to  be  controlled  in  order  to  make  AROD  fly.  These 
subsystems  are  : 

1.  Roll  rate 

2.  Rate  of  vertical  climb  . 

3.  Pitch  angle 

d.  ■  Yaw  angle 

Classical  design  tools  such  as  Root  Locus  diagrams  and  Bode  plots  are  not  well  suited 
for  MI. MO  system  design.  Instead,  the  usefulness  of  these  methods  is  primarily  liniited 
to  Single-Input  Single-Output  (SISO)  systems.  These  systems  are  generally  represented 
in  terms  of  their  S  domain  or  Z  domain  transfer  functions.  The  poles  and  zeros  of  these 
transfer  functions  determine  how  the  time  response  of  the  system  will  behave.  By  using 
the  graphical  and  analytic  methods  available  through  classical  design  theory,  the 
engineer  can  generally  place  the  poles  and  zeros  of  his  SISO  controller  in  such 
locations  as  to  obtain  an  acceptable  time  response  for  the  system.  The  complex 
interactions  that  typically  accompany  a  MIMO  system  can  become  impossible  to 
represent  in  terms  of  standard  transfer  functions.  Thus,  classical  design  methods  may 
become  powerless  for  some  systems.  By  developing  a  state  space  representation  of  the 
system,  however,  the  interactions  can  be  accurately  modelled.  Optimal  control  theory 
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is  founded  on  the  state  space  representation  of  control  systems  and,  therefore,  it  seems 
logical  to  pursue  this  theory  for  AROD.  The  basics  of  optimal  control  theory  are 
presented  in  the  following  chapter. 


II.  OPTIMAL  CONTROL  THEORY 


A.  FEEDBACK  CONTROL 
L  Why  Use  Feedback  ? 

Feedback  control  is  familiar  to  engineers  from  all  disciplines.  In  its  simpUst 
form,  feedback  control  is  nothing  more  than  using  the  present  condition,  or  "state",  of 
a  system  to  influence  its  condition  in  the  future. 

The  advantages  of  state  feedback  control  [Ref.  1:  p.97]  can  be  summarized  in 
four  points  ; 

1.  Assuming  that  controllability  and  obser\'ability  conditions  are  satisfied,  the 
transient  time  response  of  the  system  can  be  easily  controlled  and  adjusted. 

2.  The  sensitivity  of  the  system  to  plant  parameter  variation  is  reduced. 

3.  Rejection  of  disturbance  and  noise  signals  is  improved. 

4.  Steady  state  errors  may  be  eliminated  or  reduced. 

These  benefits  are  not  free.  The  penalty  for  using  feedback  control  may 
include  disadvantages  such  as : 

1.  System  complexity  increases  because  additional  sensors  may  be  required  to 
.  measure  the  feedback  states. 

2.  Sensors  contribute  to  an  increase  in  : 

■  aT  Cost 

b.  Size 

c.  Weight 

d.  Measurement  noise 

3.  Closed  loop  gain  is  generally  lower  than  open  loop  gain. 

Despite  these  potential  drawbacks,  feedback  systems  are  widely  used  in  all  engineering 
fields. 


2.  System  Classification 

Feedback  provides  a  system  with  the  ability  to  moniter  and  alter  its 
performance.  As  the  process  advances  in  time,  the  system  is  apprised  of  the  changes 
that  occur  in  its  states.  This  state  information  may  be  real  time  or  may  be  delayed  by 
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some  finite  interval  of  time.  In  general,  systems  may.  be  separated  into  two  categories 
according  to  the  nature  of  the  signals  they  process.  Theses  categories  are  : 

1.  Continuous  time. 

2.  Discrete  time. 

An  analog  electrical  circuit  is  one  example  of  the,  first  type  of  system.  In  this 
case,  the  voltage  and  current  signals  assume  values. over  a  continuum  nf  time.  That  is, 
given  any  two  instants  in  time;:  the  changing  values  of  these  signals  may  be  distinctly 
measured.  This  is  true  regardless  of.  how  closely  the awp  tihie  tinstants  occur.  Such 
system?  are  usually  described  dby.;tLxseries  of  difBsreniiai'requaitions.-  The  Xaplace 
transform  is  extremely  use&iL-iin  aHdwing  frequency 'domain  analysis  and  design  of 
continuous  time  systems.': 

A  microprocessor -bjised -system  is  an  example',  of  the  discrete  time  system.  The 
clocked  signals  in  this! -type  of  system  are  represented  by  a  sequence  of  numbers. 
Typically,  a  sequence  sampled  data,  resultsiafrom  measuring  an  analog  signal  at 
specific  intervals  in  time.  The  time !  betwetai  measurements  is  referred  to  as  the 
sampling  interval,  or  Ac..:  The  sampBhg  frequency,  f^,  is  simply  the  inverse  of  At.  For 
an  analog  sy^cm  with  aiEonrier  transfbcmtbandlimited  to  a  maximum  frequency, 
the  Nyquist frequency,'^ is  defined imJEtquation  2.1  [Ref.  3:  p.  138]. 


2f. 


max 


(2.1) 


A  general  rule  of  thumb  ibr'the  design  engineer  [Ref.  4:  p.  404]  is  to  sample  a  system 
such  that 


f,  i  10f„  (2.2) 

This  guideline  .for  selecting  a  sampling  frequency  is  based  on  the  following 
considerations 

1.  Most  systems  are  not  strictly  bandlimited.  The  choice  of  a  sampling  frequency 
greater  than  the  Nyquist  frequency  compensates  for  contamination  by  higher 
frequency  disturbances  [Ref.  5:  p.  30]. 

2.  The  sampling  frequency  should  be  fast  enough  to  avoid  aliasing.  This 
distortion  is  generated  during  the  convolution  reconstruction  of  a  time  signal 
from  its  Fourier  transform  [Ref.  3:  pp.  135-137]. 
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Discrete  time  systems  are  represented  by  difTerence  equations  and  the  Z-transform  is 
the  mechanism  by  which  these  equations  are  anaKved.  More  details  of  the 
comparisons  between  continuous  time  systems  and  discrete  time  systems  will  arise 
during  subsequent  discussion. 

3.  System  Structure  with  Feedback 


—  Figure  2.1  Basic  Control  System 

The  basic  form  for  any  control  system  is  illustrated  in  Figure  2.1.  It  is  the 
responsibility  of  the  design  engineer  to  determine  any  or  all  of  the  items  designated  in 
this  schematic.  In  this  section,  emphasis  is  placed  on  the  feedback  gain  "black  box" 
shown  in  Figure  2.1.  The  two  methods  most  commonly  used  in  practice  "to  determine 
feedback  gains  arc  pole  placement  techniques  and  optimal  control  techniques.  The 
element  of  trial  and  error  is  inherent  in  both  methods. 

Pole  placement  techniques  include  analytical  methods  as  well  as  the  frequency 
domain  methods  previously  mentioned.  These  methods  are  best  suited  to  low-order, 
linear,  time-invariant,  SISO  systems.  Although  these  methods  arc  extremely  useful  for 
certain  problems,  no  detailed  explanation  of  these  classical  techniques  is  included  in 
this  thesis. 


Optimal  control  theory  provides  an  alternative  to  classical  pole  placement 
techniques.  A  primary  advantage  of  optimal  control  methods  is  that  feedback  gains 
can  be  computed  for  a  much  broader  range  of  control  problems.  Specifically,  optimal 
control  provides  solutions  for  high  order,  non-linear,  time  varying,  MI  MO  systems. 
Such  systems  are  intractable  with  classical  methods.  In  addition,  optimal  control 
affords  the  designer  the  option  to  specify  a  performance  criteria  which  is  not  linked  to 
such  standard  time  domain  criteria  as  rise  time,  percent  overshoot,  and  settling  time. 
For  instance,  using  optimal  control  theory,  the  design  engineer  may  compute  feedback 
gains  which  result  in  a  system  that  responds  in  minimum  time  to  a  given  command 
input.  Selection  of  a  different  performance  criteria  might  result  in  a  system  that 
responds  with  minimum  energy  expendiature,  minimum  fuel,  or  minimum  deviation  from 
the  reference  command.  Sound  engineering  judgement  and  a  thorough  understanding 
of  the  system  dynamics  are  prerequisites  for  effective  application  of  optimal  control 
theory.  No  guarantee  is  made  that  the  feedback  gains  obtained  by  optimal  control 
theory  will  result  in  acceptable  system  response.  The  designer  should  evaluate  the 
system  response  and  modify  his  performance  criteria  in  order  to  achieve  the  desired 
output. 

B.  SYSTEM  DEFINITION 
I.  Continuous  Time  Systems 

The  foundation  of  a  successful  control  system  is  an  accurate  model  of  the 
plant  which  is  being  controlled.  The  state  space  representation  of  a  general  n^^  order 
continuous  time  system  is  described  by  the  following  matrix  state  equations 


x(t)  -  A(t)  x(t)  +  B(t)  u(t) 

(2.3) 

y(t)  -C(t)x(t)  +  D(t)u(t)  ■ 

(2.4) 

e(t)  -  x(t)  -  r(t) 

(2.5) 

u(t)  -  F(t)  {x(t)  -  r(t)} 

(2.6) 

where  the  definitions  in  Table  I  apply  to  a  system  with  i  control  inputs 
and  m  measurable  outputs. 
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TABLE  I 

STATE  SPACE  DEFINITIONS  FOR  CONTINUOUS  TIME  SYSTEMS 


Term 

Dimension 

Definition 

x(t) 

(n  X  1) 

State  vector 

u(t) 

(e  X  1) 

Control  input  vector  (0  <  C  2S  n) 

..  y(t) 

(m  X  1) 

Output  vector  (0  <  m  i  n) 

r(t) 

(m  X  j) 

Command  input  vector 

e(t) 

(m  X  1) 

Error  vector 

A(t) 

(n  X  n) 

Plant  matrix 

B(t) 

(n  X  e) 

Control  distribution  matrix 

C(t) 

(m  X  n) 

Output  distribution  matrix 

D(t) 

(m  X  f) 

Feedforward  control  gain  matrix 

F(t) 

(C  X  m) 

State  feedback  gain  vector 

In  this  thesis,  a  linear  time  invariant  system  will  be  assumed.  This  allows  the 
time  dependency  of  the  process  matrices,  ^(t)  and  B{t),  and  the  measurement  matrices, 
C(t)  and  D(t),  to  be  eliminated.  Because  optimal  control  theory  requires  that  all  n 
states  be  available  for  feedback,  the  output  distribution  matrix,  C,  is  set  equal  to  the 
identity  matrix,  I.  This  indicates  that  m  -  n  and  that  the  state  vector  is  completely 
observable.  In  addition,  the  feedforward  control  gain  matrix,  D,  is  assumed  to  be 
equal  to  the  zero  matrix.  These  assumptions  lead  to  the  following  simplified  state 
equations : 


x(t)  -  A  x(t)  +  B  u(t) 

(2.7) 

y(t)  -  x(t) 

(2.8) 

e(t)  -  x(t)  -  r(i) 

(2.9) 

u(t)  «  F(t)  e(t) 

(2.10) 
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Figure  2.2  Continuous  Time  System 

The  realization  of  such  a  system  is  schematically  illustrated  in  Figure  2.2. 

2.  Discrete  Time  Systems 

Optimal  control  theory  is  applicable  to  the  continuous  time  system  presented 
in  the  preceeding  section.  The  remainder  of  this  thesis,  however,  will  focus  on  the 
application  of  optimal  control  theory  to  sampled  data  systems.  The  motivation  behind 
this  effort  is  to  develop  an  interactive,  user-friendly  software  package  that  can  be 
implemented  on  a  microcomputer.  The  digital  nature  of  sampled  data  systems  make 
them  ideal  for  analysis  and  design  with  these  high  speed  computers.  The  theory  for 
optimal  control  of  discrete  time  systems  is  w'ell  developed  and  closely  follows  the 
development  for  continuous  time  systems  (Ref.  6). 

As  was  noted  earlier,  many  digital  systems  are  the  result  of  periodic  sampling 
I  of  analog  systems.  This  fact  makes  it  necessary  to  mathematically  connect  the  two 

!  types  of  systems.  For  an  analog  signal  that  is  sampled  at  the  frequency,  f ,  a  discrete 

28 

! 

I 

I 


signal  value  is  measured  evcr>' t  «  kAt  seconds.  In  this  notation,  k  is  an  integer  time 
index  in  the  range  0£  kiS  (N*l)  where  N  represents  the  last  sample  period  of  interest. 
Letting  the  sample  period  be  denoted  as 

At  -  T  (2.11) 

and  substituting  a  discrete  approximation  for  the  derivative  in  Equation  2.7, 


x(kT) 


x((k+l)T)-  x(kT) 
T 


(2.12) 


yields  the  discrete  state  equation  : 

x((k+ 1)T)  -  (I  +  AT)x(kT)  +  TBu(kT)  (2.13) 


The  analytic  solution  for  the  discrete  problem  is  given  by  : 


x(k*i-l) 

-  0x(k)  ru(k) 

(2.14) 

y(k) 

-  x{k) 

(2.15) 

e(k) 

-  x(k)  -  r(k) 

(2.16) 

u(k) 

-  F(k)e(k) 

(2.17) 

where  <I>  and  F  arc  defined  as  : 

O  -  e^"^  (2.18) 


r  -  JT  dt  B 
•'o 


(2.19) 
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The  vectors  and  matrices  which  describe  an  n*  order  discrete  time  system 
with  t  control  inputs  and  m  measured  outputs  arc  summarized  in  Table  2.  A  graphical 
realization  of  such  a  system  is  illustrated  in  Figure  2.3. 


TABLE  2 

STATE  SPACE  DEFINITIONS  FOR  DISCRETE  TIME  SYSTEMS 


*  Term 

Dimension 

Definition 

x(k) 

(n  X  1) 

State  vector 

u{k) 

it  X  1) 

Control  input  vector  (0  <  J  S  n) 

y(k) 

(m  X  i) 

Output  vector  (0  <  m  ^  n) 

r(k) 

(m  X  1) 

Command  input  vector 

e(k) 

(m  X  1) 

Error  vector 

<l»(k) 

(n  X  n) 

State  transition  matrix 

r(k) 

(n  X  e) 

Discrete  Control  distribution  matrix 

C(k) 

(m  X  n) 

Output  distribution  matrix 

D(k) 

(m  X  C) 

Feedforward  control  gain  matrix 

F(k) 

(C  X  m) 

State  feedback  gain  vector 

on  a 


The  computation  of  the  discrete  process  matrices,  <I>  and  F,  from  the 
continuous  process  matrices,  A  and  B,  is  readily  accomplished  [Ref.  5:  p.  37] 
digital  computer  as  follows. 

Define  an  auxiliary  matrix,  11  as 


n 


I  jT  gAt 

■'o 

AT^  A^T^ 
IT  + - +  - + 


2! 


3! 


(2.20) 


AiT(i+  1) 

+  .  4- 

(i+l)! 


The  terms  in  this  Taylor  series  expansion  are  computed  until  the  result  is  within  a 
specified  degree  of  accuracy.  It  behooves  the  programmer  to  set  a  very  small  tolerance 
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Figure  2.3  Time  Invariant  Discrete  Time  System 

on  the  difference  between  successive  terms  in  the  expansion  since  this  calculation  is  the 
critical  link  between  the  A  and  B  matrices  of  the  continuous  time  system  arid 
the  <I>  and  F  matrices  of  the  discrete  time  system.  Note  that  this  calculation  need 
only  be  done  once  for  a  given  system  with  a  fixed  sample  interval,  T.  The  link  is 
completed  by  using  Equations  2.21  and  2.22. 

O  -  I  +  An  (2.21) 


r  -  OB  (2.22) 

The  subroutine  PHIDEL  listed  in  Appendix  C  performs  the  calculations  required  to 
solve  Equations  2.20,  2.21,  and  2.22.  A  tolerance  of  10*^  is  used  for  the  allowable 
error  between  the  last  term  used  from  the  Taylor  e.\pansion  and  the  first  term  not  used. 


3.  Constraints 

The  system  is  now  defined  in  terms  of  the  O  and  F  state  space  matrices.  The 
next  step  in  the  design  process  is  to  identify  any  constraints  under  which  the  system 
will  operate.  These  constraints  are  as  unique  to  the  problem  at  hand  as  is  the  system 
itself.  Therefore,  no  detailed  discussion  on  constraints  is  appropriate  without  first 
defining  a  specific  problem.  This  is  done  in  Chapter  III. 

C.  THE  PERFORMANCE  MEASURE 

1. *  Quadratic  Cost  Function  "■ 

The  central  theme  in  discrete  optimal  control  theory  is  minimization  of  a  cost 
function,  J,  defined  in  Equation  2.23. 

J  -  e‘(N)He(N)  +  ^  [e‘(k)Qe(k)  +  u‘(k)Ru(k)]  (2.23) 

where 

J  -  Scaler  cost  of  operatina  the  system  over 
the  time  mterval  0:S  fs  (N-1) 

e<N)  *  Slate  vector  at  the  terminal  time 

e(k)  “  State  veaor  at  intermediate  discrete  times 

u(k)  -  Control  vector  at  intermediate  discrete  times 

_  H  *•  Terminal  state  weighting  matrix 

Q  -  State  trajectory  weighting  matrix 

R  *■  Scaler  control  weighting  matrix 

N  •  Time  index  at  terminal  lime 

t  -  Matrix  transpose  operator 

2.  Regulators  and  Trackers 

It  is  imperative  to  note  here  that  the  error  state  vector,  e(k),  in  Equation  2.23 
may  not  be  the  same  as  the  system  state  vector,  x(k),  that  is  presented  in  the  previous 
section.  The  e  state  vector  is  usually  formulated  in  one  of  two  ways  : 

1.  If  e(k)  -  x(k),  then  the  problem  is  a  'regulator'  problem.  The  objective  is  to 
drive  the  system  states  to  the  origin  during  the  time  interval  (1,N). 

2.  If  e(k)  -  x(k)  -  r(k),  then  the  problem  is  a  'tracking'  problem.  The  objective  is 
to  drive  the  system  states  to  have  minimum  deviation  from  the  conunand  input, 
r(k),  during  the  time  interval  (I, N). 
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In  order  to  extend  the  regulator  problem  to  the  tracking  problem,  the  command  input 
vector,  r(k),  must  contain  the  same  eigenvalues  and  structure  as  the  x(k)  state  vector. 
It  is  possible,  in  many  problems,  to  structure  an  approximate  r(k)  vector  so  that  the 
use  of  the  regulator  solution  is  allowed.  In  the  case  of  the  regulator,  the  control  input 
signal  is  generated  as  shown  in  Equation  2.24. 

u(k)  -  F(k)  x(k)  (2.24) 


In  the  case  of  the  tracking  problem,  the  control  signal  is  generated  from  the  error 
signal  as  shown  in  Equation  2.25. 

u(k)  -  F(k){x(k).r(k)}  (2.25) 


The  comparison  between  these  two  types  of  systems  is  demonstrated  in  their  block 
diagram  structures  as  shown  in  Figure  2.4. 

3.  Performance  Weighting  Factors 

The  H,  Q,  and  R  weighting  matrices  are  the  parameters  by  which  the  design 
engineer  shapes  the  solution  of  an  optimal  control  problem  to  best  suit  the  problem. 
There  are  no  magic  formulas  which  instruct  the  designer  on  how  to  choose  the  values 
of  these  parameters.  Intuition,  experience,  and  patience  are  the  keys  to  successful 
design.  It  is  in  selecting  values  for  these  performance  weighting  factors  that  the 
process  of  trial  and  error  enters  the  design  process.  There  are,  however,  some 
restrictions  and  general  guidelines  that  partially  direct  the  efforts  of  a  design  engineer. 

First  consider  the  restrictions.  Both  of  the  state  weighting  matrices,  H  and  Q, 
must  satisfy  all  of  the  criteria  below  [Ref.  7;  p.  753]. 

1.  Dimension  is  (n  x  n) 

2.  Real 

3.  Symmetric 

4.  Positive  semidefmite 

In  addition,  the  designer  should  never  allow  both  H  and  Q  to  be  equal  to  the  zero 
matrix  at  the  same  time.  The  resulting  cost  function  would  be 
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Figure  2.4  Comparison  Between  Regulator  and  Tracker  System  Structure 


(2.26) 


N-l 

J  ■*  Z  «*(k)Ru(k) 

k»0 

It  is  simple  to  see  that  the  cost  function  will  be  minimized  by  setting  u(k)  equal  to  zero 
[Ref.  7:  p  755].  This  would  be  disasterous  since  there  would  be  no  command  signal 
available  to  drive  the  system  states  toward  the  desired  state.  It  is  permissible  to  set 
either  H  or  Q  equal  to  the  zero  matrix  provided  that  they  are  never  both  zero 
simultaneously.  The  {t  x  C)  control  weighting  matrix,  R,  must  be  positive  defmite  in 
order  to  assure  the  existence  of  a  finite  control  [R.ef.  7:  p.  754]. 

Although  there  is  no  requirement  for  H  and  Q  to  be  diagonal  matrices,  the 
usual  practice  is  to  select  non-negative  values  for  their  diagonal  elements  and  to  set  all 
ofT-diagonal  elements  to  zero.  This  technique  eliminates  all  cross  product  terms  in  the 
cost  function.  For  example,  consider  a  second  order  system  containing  states  Xj  and 
X2.  If  diagonal  matrices  are  selected  for  H  and  Q,  then  only  the  (x^)^  and  (Xj)^  terms 
will  be  weighted  in  the  cost  function.  There  will  be  no  consideration  given  to  the  XjX2 
or  X2Xj  cross  product  terms. 

Assuming  that  neither  H  nor  Q  is  the  zero  matrix,  it  is  suggested  that  the 
elements  of  these  weighting  matrices  be  selected  so  as  to  normalize  the  values  of  the 
states  which  they  multiply  [Ref.  6:  p.  32].  For  instance,  suppose  that  Xj  represents  the 
RPM  of  AROD'sqjropeller  and  X2  represents  the  angular  displacement  in  degrees  of 
the  throttle  servo.  State  Xj  is  expected  to  have  a  nominal  value  of  7200  while  state  X2 
may  have  a  nominal  value  of  only  10.  Assume  that  element  q^j  which  multiplies  (Xj)^ 
and  element  q22  which  multiplies  1X2)^  are  both  set  equal  to  one  in  an  attempt  to 
weight  the  two  states  equally.  In  terms  of  the  cost  function,  the  RPM  will  be  weighted 
much  heavier  (approximately  720^  times  heavier  !!)  than  the  throttle  servo  position 
angle.  An  appropriate  solution  is  to  set  qjj  (1/720)^  if  ^22'^  Thus,  each  signal 
is  given  equal  consideration  in  the  cost  function. 

Notice  in  Equation  2.23  that  the  terminal  cost  term  is  not  included  inside  the 
Y,  operator.  This  means  that  the  H  term  is  only  counted  one  time  and  therefore  can 
contribute  only  once  to  the  overall  cost.  The  state  trajectory  term  and  the  control 
term,  however,  are  counted  N  times.  If  no  compensation  is  made  for  this  disparity,  the 
cost  of  an  error  in  the  terminal  state  is  likely  to  not  have  any  impact  on  the  control 
solution.  It  seems  that  an  additional  scaling  is  required  on  the  weighting  elements.  If 
the  summation  in  Equation  2.23  is  to  be  made  over  500  discrete  time  intervals,  for 
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instance,  the  normalized  elements  of  H  should  be  multiplied  by  500  in  order  to  weight 
the  terminal  states  on  the  same  scale  as  the  state  trajectoiy  and  control  eflbrt  terms. 
Alternately,  the  elements  of  Q  and  R  could  be  divided  by  500.  It  is  the  relative 
magnitudes,  not  the  absolute  magnitudes,  of  these  weighting  factors  which  shape  the 
control  solution. 

4.  Speciflc  Types  of  Problems 

Optimal  control  theory  can  be  used  to  solve  any  of  the  following  types  of 
problems : 

1.  Minimal  time  ~ 

2.  Minimal  control  effort 

a.  Minimum  fuel 

b.  Minimum  energy 

3.  Minimal  error  from  a  reference 

a.  Regulator 

b.  Tracking 

c.  Terminal  state  control 

Each  of  these  problem  types  requires  minimization  of  a  unique  cost  function  in  order 
to  generate  the  appropriate  control  solution  [Ref.  6:  pp,  30-34].  The  cost  functions 
which  are  to  be  explored  in  this  thesis  are  listed  in  Table  3. 


TABLE  3 

TYPICAL  COST  FUNCTIONS 

Goal 

Cost  Function 

Additional  Explanation 

Regulator 

Jl  -  I  x^k)Qx(k) 

Tracker 

“  I  e'(k)Qe(k) 

e(k)  -  x(k).r(k) 

Terminal  Control 

J3  -  eVN)He(N) 

e(N)  -  x(N)  -  r(N) 

Minimum  Energy 

J4  •  X  u^k)Ru(k)  +  J3 

Minimum  Fuel 

J5  -  au(k)|. 

All  Y,  af*  performed  over  the  interval  0  k  ^  (N-1). 
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Of  course,  the  cost  function  in  Equation  2.23  is  a  general  form  which  embodies  all  of 
the  specific  cost  functions  (except  for  minimum  fuel)  contained  in  Table  3.  Proper 
selection  of  H,  Q,  and  R  will  allow  calculation  of  feedback  control  gains  which  cause 
the  system  to  meet  the  specified  performance  goal. 

D.  CALCULATION  OF  OPTIMAL  FEEDBACK  GAINS 

The  method  of  dynamic  programming  is  the  workhorse  which  permits  the 
calculation  of  optimal  feedback  control  gains.  Developed  in  the  late  1950's  by  R.E. 
Bellman,,  this  design  tool  provides  a  closed  form  solution  for  the  minimization  j)f  the 
quadratic  cost  function  for  a  discrete  time  linear  system  [Ref.  6:  p.  84]. 

The  procedure  involved  in  calculating  optimal  control  gains  is  unique  in  that  the 
computation  begins  with  the  final,  or  N*,  discrete  time  interval  and  progresses 
backwards  in  time  to  the  previous  interval  of  the  system  process.  This  procedure  in 
'negative  time'  is  possible  because  the  calculation  of  the  gains  do  not  require  any 
information  about  the  state  vector,  x(k).  The  sequence  of  calculations  continues  in 
negative  time  until  a  separate  gain  matrix  is  determined  for  every  discrete  sample 
period  in  the  time  interval  (0,N).  The  resulting  time-dependent  gain  trajectory  is 
usually  stored  in  memory  so  that  the  control  signal,  u(k),  may  be  computed. 

An  interesting  and  very  useful  property  of  the  gain  trajectory  is  its  tendency  to 
approach  a  constant  valued  matrix,  under  certain  conditions  [Ref.  4:  p.  354].  This 
constant  matrix  is  referred  to  as  the  steady  state  feedback  gain  matrix.  The  conditions 
necessa^  for  F^^  to  exist  are  : 

1.  The  system  is  controllable. 

2.  the  and  F  Matrices  are  time  invariant. 

3.  The  H  terminal  state  weighting  matrix  is  equal  to  the  zero  matrix. 

4.  The  Q  trajectory  state  weighting  matrix  is  constant. 

5.  The  R  control  weighting  matrix  is  constant. 

6.  The  number  of  stages,  N,  of  the  process  is  large. 

It  is  possible  for  the  feedback  gain  matrix  to  approach  F^^  without  satisfying  the  first 
three  conditions  above.  When  all  conditions  are  satisfied,  however,  a  steady  state  gain 
solution  is  guaranteed  provided  that  N  is  large  enough.  Just  how  large  M  must  be  in 
order  to  allow  the  gain  trajectory  to  reach  steady  state  is  determined  by  the  slowest 
time  constant  in  the  solution  of  the  discrete  matrix  Riccati  equation  [Ref.  5:  p.  259].  In 
practice,  trial  and  enor  is  the  the  most  expedient  method  to  determine  how  large  N 
must  be  in  order  to  ensure  a  steady  state  gain  matrix,  F^^. 
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A  series  of  three  recursive  equations  (Ref.  6:  p.  83]  is  used  to  calculate  the  gain 
trajector>’,  F(k).  It  is  convenient  to  introduce  a  negative  time  discrete  index,  K,  which 
is  defined  as  follows 


K  -  N  -  k  (2.27) 

Since  k  varies  from  (0,N-1)  as  fonvard  time  progresses,  K  varies  from  (1,N)  as  negative 
time  progresses.  Equation  2.28  is  the  solution  for  the  transpose  of  the  optimal 
feedback  gain  vector  at  each  discrete  time  step.  This  equation  is  the  solution  to  the 
well  known  discrete  matrix  Riccati  equation.  Equations  2.29  and  2.30  are  auxiliary 
equations  required  to  complete  the  calculations.  I  he  recursive  matrix  equations  are  : 


F(K)  -  -  {r‘  P(K-l)  r  +  R}-*  {F^  P(K-l)  a>}  (2.28) 

*P(K)  -  O  +  r  F(K)  (2.29) 

P(K)  -  T'(K)  P(K-I)  T(K)  +  Q  +  F^K)R  F(K)  (2.30) 

with  'negative  time'  initial  conditions 

F*(0)  «  0  (2.31) 

T(())  -  0  (2.32) 

P(0)  -  H  (2.33) 


While  somewhat  difficult  at  first  glance,  Equations  2.28,  2.29,  and  2.30  are  ideal  for 
iterative  solution  by  a  digital  computer.  These  equations  arc  solved  in  the  main 
OPTCON  program  listed  in  Appendix  B.  The  reader  who  is  not  familiar  with 
OPTCO.N  is  encouraged  to  review  the  discussion  of  this  program  in  Appendix  A  prior 
to  continuing  with  Chapter  III. 


III.  CONTROL  SYSTEM  DESIGN  FOR  AROD 


A.  OVERVIEW 

The  purpose  of  this  chapter  is  to  use  optimal  control  theory  to  design  an 
automatic  flight  control  system  for  the  U.S.  Marine  Corps'  remotely  piloted  AROD. 
In  their  initial  form,  the  equations  of  motion  which  describe  AROD's  dynamic 
behavior  are  extremely  nonlinear  and  present  a  formidable  challenge  to  the  control 
system  designer.  For  this  reason,  the  equations  are  first  linearized  about  a  steady  state 
hover  condition.  The  restrictions  and  assumptions  under  which  the  linearized 
equations  of  motion  are  developed  are  summarized  below: 

1.  The  mass  of  AROD  is  constant  with  time  [Ref.  8:  p.  11], 

2.  The  propeller  angular  velocity  is  constant. 

3.  Perturbations  from  steady  state  hover  are  small.  This  restriction  requires  that 
AROD  pitch,  roll,  and  yaw  angular  displacements  be  limited  to  less 
than  15«  (Jt/12  *  0.2618  radians)  (Ref.  8:  p.  41]. 

4.  Steady  state  pitch  and  yaw  rates  are  zero. 

5.  Initial  side  velocities  are  zero. 

6.  Initial  bank  angles  are  zero. 

7.  Initial  angular  velocities  are  zero  [Ref.  8:  p.  45). 

In-order  to  elucidate  the  equations  of  motion,  Figure  3.1  is  provided  for 
reference.  The  axis  system  in  Figure  3.1  is  known  as  a  body-fixed  coordinate  system. 
The  body-fixed  axis  system  is  thought  of  as  being  rigidly  attached  to  the  vehicle  so  that 
any  change  in  the  orientation  of  AROD  with  respect  to  the  earth-fixed  axis  system 
(X',  Y',  Z')  results  in  a  corresponding  change  in  the  orientation  of  the  body-fixed  axis 
system  (X,  Y,  Z)  with  respect  to  (X',  Y',  Z').  The  angles  <p,  0,  and  y,  commonly 
known  as  the  Euler  angles  [Ref.  8:  p.  25],  are  measures  of  the  roll,  pitch,  and  yaw 
angles  respectively  between  the  body-fixed  (X,  Y,  Z)  coordinate  system  and  the  earth- 
fixed  <X',  Y',  Z')  coordinate  system.  The  angular  velocities  p,  q,  and  r  in  Figure  3.1 
correspond  to  the  roll,  pitch,  and  yaw  rates  respectively. 

The  automatic  flight  control  system  is  logically  separated  into  three  subsystems 
according  to  the  simplified  equations  of  motion.  The  three  control  subsystems  which 
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Figure  3.1  AROD  Body- Fixed  Coordinate  System 

are  hereaf^r  designed  are  : 

1.  Roll  rate  controller. 

2.  Altitude  rate  controller. 

3.  Pitch  angle  and  yaw  angle  controller. 

Each  of  these  controllers  is  designed  independently  from  the  other  subsystems. 
Therefore,  any  cross-coupling  which  may  occur  across  the  subsystem  boundaries  is  not 
accounted  for.  Each  of  the  following  sections  is  devoted  to  the  design  of  a  controller 
for  one  specific  AROD  subsystem.  A  detailed  design  is  first  presented  in  Section  B  for 
the  roll  rate  controller  in  order  to  demonstrate  the  iterative  nature  of  the  design 
process.  Section  C  presents  the  results  for  the  altitude  rate  controller.  In  the  interest 
of  brevity,  only  the  initial  trial  run  and  the  final  solution  for  the  altitude  rate  controller 
are  presented.  The  couphd  dynantics  ol  A  ROD'S  pitch  and  yaw  is  examined  in  greater 
detail  in  Section  D. 


B.  AROD  ROLL  RATE  CONTROLLER 

1.  The  Roll  System 

The  purpose  for  roll  rate  control  of  the  AROD  is  to  allow  the  remote  pilot  to 
command  a  desired  rotation  velocity  about  the  vehicle's  longitudinal,  or  x,  axis.  Such 
movement  allows  the  camera  aboard  the  vehicle  either  to  slowly  scan  a  selected  ground 
sector  or  to  terminate  the  rolling  motion  so  that  a  target  of  interest  can  be  further 
studied.  The  nature  of  remote  sensing  requires  that  the  vehicle  respond  rapidly  to  a 
roll  rate  command.  When  the  remote  pilot  locates  a  ground  target,  he  needs  to  be  able 
to  swiftly  hring  the  vehicle  to  a  zero  roll  rate  condition  with  negligible  overshoot?  Such 
movement  is  commanded  through  a  twistable  handgrip  control  located  on  the  pilot's 
console.  It  is  assumed  that  this  roll  rate  command,  p^,  is  limited  to  a  step  input 
of  1  radian/second  (57.5«/second).  Although  no  time  response  criteria  are  specified  by 
the  Marine  Corps,  it  is  assumed  for  the  purpose  of  this  work  that  the  following  design 
specifications  for  roll  rate  are  required  : 

1.  Zero  steady  state  error  for  a  step  input  is  required. 

2.  The  two  percent  settling  time,  t2%’  ^  second. 

3.  No  overshoot  is  allowed. 

The  simplified  equation  of  motion  which  describes  AROD's  roll  rate 
subsystem  is  given  in  Equation  3.1. 

% 

P-L.S.  (3.1) 

The  aileron  servo  dynamics  are  modelled  in  Equation  3.2  as  a  second  order  plant  with 
a  natural  frequency,  oi,  of  2  Hz  and  a  damping  coefficient,  of  0.707. 

•  •  '  • 

dj,  »  *2IJ(adjj  -  +  oj^Ug  (3.2) 

The  definitions  in  Table  4  apply  to  Equations  3.1  and  3.2.  In  order  to  apply  the  theory 
of  optimal  control,  a  suitable  state  space  representation  of  the  roll  rate  system  must 
first  be  developed.  Figure  3.2  presents  the  state  space  signal  flow  graph  selected  for 
this  subsystem. 


TABLE  4 

VARIABLE  DEFINITIONS  FOR  AROD  ROLL  RATE  EQUATIONS  OF 

MOTION 


Variable 

Definition 

Value 

Units 

P 

Vehicle  Roll  Rate 

TBD 

radians/ second 

Aiierop  Effectiveness 
Coefficient 

-21.29 

seconds'^ 

«. 

• 

AileroA  Servo 

Deflection  Angle 

^  30o 

radians 

Aileron  Servo 

Deflection  Velocity 

i  50»/sec 

radians/second 

a 

; 

Ailerop  Servo 

Dampmg  CoefHcient 

0.707 

unitless 

CO 

Aileron  Servo 

Natural  Frequency 

12.57 

radians/second 

“a 

Contfol  Input 
to  Aileron  Servo 

TBD 

volts 

By  designating  the  output  of  each  integrator  in  Figure  3.2  to  be  a  state,  the  following 
third  order  state  equations  are  derived  : 


(3.3) 


0 

-21.29 

ol 

0 

0 

1 

0 

-157.91 

-17.77| 

(3.4) 


u^  »  F  {x  -  r} 


(3.5) 


Assuming  that  a  unit  step  roll  rate  is  commanded,  the  command  input  vector  becomes 


m 

' 

•  ■ 

Pc 

a 

r  * 

|ac 

ai 

0 

®ac 

0 

(3.6) 
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where  the  subscipt  c  indicates  a  command  input  variable.  Now  that  the  roll  rate 
system  and  its  constraints  have  been  identified,  the  next  step  is  to  use  the  OPTCON 
program  to  design  a  suitable  roll  rate  controller. 

2.  Roll  Rate  Controller  Design 
a.  Choosing  a  Sampling  Frequency 

In  order  to  determine  an  appropriate  sampling  rate  for  the  roll  rate  system 
given  in  Equations  3.3  through  3.6,  the  bandwidth  of  the  open  loop  system  is  first 
determined.  The  open  loop  transfer  function  for  the  roll  rate  system  is  given  in 
Equation  3.7.  " 


P(s)  ^  (157.91)  (21.29) 

U^(s)  ’  “  s^  +  17.77  s  +  157.91  ’ 

Tht  open  loop  Bode  diagram  for  this  transfer  function  is  shown  in  Figure  3.3.  The 
negative  3  dB  bandwidth  of  the  magnitude  curve  is  approximately  12 
radians  per  second  or  2  Hz.  Using  the  criterion  discussed  in  Chapter  2,  The  Nyquist 
sampling  rate,  f^,  is  4  Hz.  In  order  to  avoid  aliasing  effects  and  to  ensure  that  the 
sampled  system  is  a  reasonable  representation  of  the  continuous  time  system,  it  is 
decided  to  employ  a  sampling  frequency  that  is  at  least  five  times  greater  than  f^^.  A 
comparison  of  the  effect  of  using  various  sampling  frequencies  is  given  in  Table  5.  The 
cost  function  which  is  used  to  obtain  the  optimal  gains  for  this  comparison  is  included 
in  Table.  5  and  is  hereafter  referred  to  as  the  baseline  cost  function.  The  colunnn 
labelled  "Number  of  Stages  Required"  in  Table  5  refers  to  the  number  of  discrete  time 
stages  that  must  elapse  before  the  optimal  gain  vector  reaches  its  steady  state  value, 
F,,.  Notice  that  the  magnitude  of  the  steady  state  gain  vector  is  related  to  the 
sampling  frequency.  In  general,  a  faster  sampling  rate  yields  larger  feedback  gains. 
Also  notice  that  the  sampling  rate  affects  the  amount  of  real  time  that  is  required  for 
the  gains  to  reach  steady  state.  For  example,  when  f,  is  20  Hz,  it  requires  1.60  seconds 
for  Fjj  to  be  achieved.  When  f^  is  40  Hz,  however,  a  total  time  of  2.23  seconds  elapses 
before  F^,  is  achieved.  These  considerations  are  important  if  the  control  gains  are  to 
be  dynamically  implemented.  In  the  case  of  this  design,  the  control  implementation  is 
limited  to  steady  state  values  of  the  optimal  gains.  The  unit  step  time  response 
obtained  by  using  steady  state  optimal  feedback  gains  is  observed  to  be  nearly  identical 
for  all  three  of  the  sampling  rates  listed  in  Table  5.  Specifically,  the  roll  rate  state,  x^ 
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Figure  3.3  Open  Loop  Bode  Diagram  For  Roll  Rate  System 
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TABLE  5 

EFFECT  OF  SAMPLING  FREQUENCY 
ON  ROLL  RATE  SYSTEM 


Baseline  Cost  Function 

■“■[ill] 


time  response  exhibits  an  overshoot  of  approximately  3.7%  and  a  2%  settling  time  of 
1.3  seconds.  See  Figure  3.4.  This  implies  that  any  of  the  three  sampling  rates 
examined  is  acceptable.  Because  it  is  generally  considered  good  practice  to  implement 
small  feedback  gains  when  possible,  it  is  decided  to  employ  a  sampling  frequency  of  20 
Hz  for  the  remainder  of  the  roll  rate  controller  design. 
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b>  Methodology 

At  this  point,  four  different  groups  of  cost  functions  are  examined.  For  a 
given  cost  function  group,  the  H  and  Q  matrices  are  held  constant  while  the  control 
weighting  factor,  R,  is  varied  within  the  range  (0.01,  100).  Approximately  15  runs  are 
made  for  each  cost  function  group  with  a  different  value  of  R  inserted  into  the  cost 
function  for  each  run.  After  the  steady  state  gains  are  determined  for  a  given  run,  they 
are  implemented  into  the  control  equation  and  a  time  response  for  the  roll  rate  state, 
x^,  is  obtained.  The  percent  overshoot  and  2%  settling  time  are  recorded  for  each  Xj 
time  response.  A  summary  of  this  information  is  presented  for  the  four  cost  function 
groups  in  Tables  6,  7,  8,  and  9.  Following  each  of  these  four  tables,  there  appears  a 
graph  of  two  time  response  parameters,  percent  overshoot  and  settling  time,  versus  the 
value  of  the  control  weighting  factor,  R. 

It  is  hardly  necessary  to  include  a  time  response  graph  for  all  of  the  59 
total  runs  examined.  It  is  instructive,  however,  to  compare  the  time  responses  for  a 
selected  set  of  cost  functions.  Three  runs  in  each  of  the  parameter  summary  tables. 
Tables  6  through  9,  are  flagged  with  asterisks.  These  flags  indicate  that  the  roll  rate 
time  response  graph  for  that  run  is  included  subsequent  to  the  applicable  summary 
table.  Keep  in  mind  that  the  criteria  for  the  roll  rate  step  response  is  specified  to  be 
such  that  there  is  no  overshoot  and  the  2%  settling  time  is  less  than  one  second.  A 
discussion  of  the  results  from  the  four  cost  function  groups  follows  the  last  figure  in 
this  series  of  tables  and  graphs. 


TABLE  6 

ROLL  RATE  PARAMETERS  FOR  GROUP  1 


Q  ■■ 


Sampling  Interval  T 

■  0.05  seconds 

Run 

Control 

We^ht 

fl 

Steady  State 

Gams 

Percent 

Overshoot 

Settling 

Time 

(sec) 

1 

0.01 

0.1727 

-0.2944 

-0.0892 

4.03 

1.31 

2 

0.03 

0.1726 

-0.2942 

-0.0891 

4.02 

1.31 

3 

0.05 

0.1725 

-0.2940 

-0.0890 

4.f2 

1.31 

4  * 

0.10 

0,1724 

-0.2936 

-0.0890 

4.00 

1.31 

5 

0.30 

0.1717 

-0.2920 

-0.0884 

3.95 

1.31 

6 

0.50 

0.1711 

-0.2904 

-0.0879 

3.89 

1.31 

-  T 

1.00 

0.1696 

-0.2866 

-0.0866 

3.74 

1.30 

8 

3.00 

0.1638 

-0.2728 

-0.0820 

3.17 

1.30 

9 

5.00 

0.1587 

-0.2610 

-0.0780 

2.68 

1.29 

10 

7.00 

0.1541 

-0.2509 

-0.0744 

2.22 

1.24 

11 

10.00 

0.1480 

-0.2380 

-0.0697 

1.65 

0.87 

12  *• 

20.00 

0.1323 

-0.2078 

-0,0582 

0.42 

1.01 

13 

30.00 

0.1209 

-0.1886 

-0;0504 

0.00 

1.17 

14 

50.00 

0.1053 

-0.1646 

-0.0403 

0.00 

1.47 

15  *•* 

100.00 

0.0835 

-0.1350 

-0.0278 

0.00 

2,05 
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Figure  3.8  Roll  Rate  Time  Response  for  Group  1  Run  Number  15 
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TABLE  7 

ROLL  RATE  PARAMETERS  FOR  GROUP  2 


« -  [1 '!  il]  ^  “I  o.f] 

Sampling  Interval  T  *  0.05  seconds 


Run 

Control 

we^ht 

fl 

Steady  State 

Gams 

h  ^3 

Percent 

Overshoot 

Sejtling 

Tune 

(sec) 

1 

0.0  i 

0.4803 

-1.2269 

-0.1073 

4.37 

0.79 

2 

0.03 

0.4786 

-1.2218 

-0.1068 

4.35 

0.79 

3 

0.05 

0.4769 

-1.2168 

-0.1063 

4.33 

0.79 

4  • 

0.10 

0.4728 

-1.2047 

-0.1052 

4.23 

0.79 

5 

0.30 

0.4576 

-1.1598 

-0.1009 

4.08 

0.79 

6 

0.50 

0.4441 

-1.1202 

-0.0971 

3.88 

0.79 

1 

1.00 

0.4159 

-1.0381 

-0.0893 

3.49 

0.79 

8 

3.00 

0.3442 

-0.8374 

-0.0700 

2.12 

0.73 

9 

5.00 

0.3024 

-0.7251 

-0.0593 

1.12 

0.57 

10 

7.00 

0.2737 

-0.6504 

-0.0522 

0.44 

0.62 

11  •* 

10.00 

0.2435 

-0.5735 

-0.0450 

0,00 

0.70 

12 

.30.00 

0.1596 

-0.3702 

-0.0268 

0.00  . 

1.16 

13 

50  00 

0.1281 

-0.2969 

-0.0206 

0.00 

1.46 

14  ••• 

100.00 

0.0937 

-0  2175 

-0.0144 

0.00 

2.00 
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Figure  3.11  Roll  Rate  Time  Response  for  Group  2  Run  Number  11 
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TABLE  8 

ROLL  RATE  PARAMETERS  FOR  GROUP  3 


«-[i 

Sampling  Interval  T  *  0.05  seconds 


Run 

Control 

we^ht 

Steady  State 

f  f 

Percent 

Overshoot 

Settling 

1 

0.01 

1.1983 

-2.6905 

-0.1332 

4.53 

0.48 

2 

0.03 

1.1622 

-2.6141 

-0.1296 

4.54 

0.49 

3 

0.05 

1.1300 

-2.5460 

•0.1264 

4.58 

0.49 

4* 

O.iO 

1.0625 

-2.4028 

-0.1196 

4.65 

0.49 

5 

0.30 

0.8917 

-2.0369 

-0.1023 

4.63 

0.51 

6 

0.50 

0.7917 

-1.8200 

-0.0919 

4.29 

0.52 

7 

1.00 

0.6497 

-1.5079 

-0.0770 

3.64 

■  0.53 

8 

1.50 

0.5689 

-1.3280 

-0.0683 

2.87 

0.54 

-t 

2.00 

0.5144 

-1,2053 

-0.0623 

2.33 

0.53 

10 

3.00 

0.4426 

-1.0425 

-0.0543 

1.28 

0.43 

11 

5.00 

0.3621 

-0.8576 

-0.0452 

0.00 

0.50 

12 

10.00 

0.2712 

-0.6458 

-0.0345 

0.00 

0.71 

13 

15.00 

0.2273 

-0.5427 

-0.0292 

0.00  ' 

0.87 

14  •• 

20.00 

0.2001 

-0.4782 

-0.0258 

0.00 

0.97 

15 

30.00 

0.1663 

-0.3986 

-0.0217 

0.00 

1.16 

16 

50,00 

0.1316 

-0.3153 

-0.0172 

0.00 

1.46 

17  *•* 

100.00 

0,0950 

-0.2277 

-0.0126 

0.00 

2.00 
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Figure  3.13  Group  3  Time  Response  Parameters 
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Figure  3. 14  Roll  Rate  Time  Response  for  Group  3  Run  Number  4 
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Figure  3.15  Roll  Rate  Time  Response  for  Group  3  Run  Number  14 
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Figure  3.16  Roll  Rate  Time  Response  for  Group  3  Run  Number  17 
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TABLE  9 

ROLL  RATE  PARAMETERS  FOR  GROUP  4 


m 


m 


Sampling 

Interval  T 

■  0.05  seconds 

Run 

Control 

We^ht 

fl 

Steady  State 

Gams 

*2  *3 

Percent 

Overshoot 

S^tling 

Tune 

(sec) 

1 

0.01 

3.1663 

-5.4007 

-0.1713 

9.57 

0.27 

2 

0.03 

2.3689 

-4.3357 

-0.1480 

8.82 

0.31 

3 

0.05 

2.0386 

-3.8586 

-0.1366 

8.67 

0.33 

4  • 

0.10 

1.6396 

-3.2461 

-0.1209 

7.37 

0.36 

.  5 

0.30 

1.1264 

-2.3821 

-0.0962 

6.22 

0.42 

6 

0.50 

0.9350 

-2.0316 

-0.0852 

5.96 

0.44 

-  J7 

1.00 

0.7182 

-1.6110 

-0,0709 

4.51 

0.48 

8 

3.00 

0.4613 

-1,0739 

-0.0506 

1.33 

0.41 

9 

5.00 

0.3718 

-0.8759 

-0.0425 

0.00 

0.49 

10 

10.00 

0.2750 

-0.6551 

-0.0329 

0.00 

0.72 

'  li 

30.00 

0.1674 

-0.4020 

-0.0210 

0.00 

1.16 

12 

50.00 

0.1320 

-0.3175 

-0.0168 

0.00 

1.46 

13 

100.00 

0.0951 

-0.2289 

-0.0123 

0.00 

2.00 
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Figure  3.17  Group  4  Time  Response  Parameters 


65 


CONTROL  WEIOIITING  FACTOR  (10 


Figure  3.18  Roll  Rate  Time  Response  for  Group  4  Run  Number  4 
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Figure  3.20  Roll  Rate  Time  Response  for  Group  4  Run  Number  13 
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c.  Resuits 


(1)  Group  i  Cost  Functions.  All  ihree  states  are  weighted  with  a  value  of 
unity  for  the  Group  1  cost  functions.  This  implies  that  the  designer  is  attempting  to 
minimize  the  error  in  all  states  with  equal  emphasis.  Pursuit  of  tliis  group  of  cost 
functions  is  made  in  order  to  determine  general  patterns  of  cause  and  effect.  For 
example,  it  is  apparent  in  Table  6  that  the  control  weighting  factor,  R,  significantly 
affects  the  magnitude  of  the  steady  state  optimal  feedback  gain  vector,  By 
increasmg  the  penalty  on  the  control  effort,  the  magnitudes  of  the  feedback  gains  are 

A  b 

reduced.  Thus,  if  the  control  system  is  physicallv  limited  to  some  maximum  value  of 
control  effort,  then  the  R  term  is  the  logical  parameter  to  adjust.  The  percent 
overshoot  and  settling  time  data  frum  Table  6  indicates  that  R  directly  affects  the  time 
response  as  well.  From  Figure  3.5,  note  that  the  value  of  R  has  negligible  effect  on  the 
time  response  parameters  foi  any  value  of  R  less  than  unity.  Refer  to  Figure  3.6  in 
tvhich  R  «  0.1.  As  the  control  weighting  factor  inn  eases  above  unity,  however,  the 
time  response  is  dramatirally  affected.  For  values  of  P,  greater  than  30,  the  time 
response  exhibits  no  overshoot  and  the  settling  time  appears  to  lengthen  without 
bovmd  as  the  system  becomes  increasingly  slow.  Refer  to  Figure  3.8  in  which 
R  -  100.  Also  notice  in  Figure  3.5  that  there  is  no  cost  function  in  Group  1  that 
yields  an  acceptable  time  response  which  satisfies  both  of  the  roll  rate  criteria.  The 
cost  function  in  Group  1  which  yields  a  timu  response  closest  to  the  design 
specifications  is  run  number  12  shown  in  Figure  3.7.  This  run  is  subsequently  used  as 
a  basIsTor  comparison  of  the  time  responses  generated  by  the  other  three  groups  of 
cost  functions. 

(2)  Group  2  Cost  Functions.  Because  acceptable  results  are  not  obtained 
from  the  Group  1  cost  fur.ciicns,  it  is  decided  to  place  increased  emphasis  on  the  error 
ill  the  x^  state  while  reducing  the  emphasis  on  the  error  in  the  x^  and  x^  states.  This 
tactic  is  allowable  because  the  maximum  absolute  values  of  the  X2  and  states  are 
significantly  less  than  the  constraints  for  the  6^^  and  6^  servo  states  listed  in  Table  4. 
Table  7  summarizes  the  data  for  Group  2.  Figure  3.9  evidences  the  relationship 
between  the  time  response  parameters  and  the  control  weighting  factor,  R.  Notice 
in  tills  figure  that  an  acceptable  time  response  is  expected  for  any  Group  2  cost 
function  in  which  10  :S  R  i  20,  Figure  3.11  shows  the  x,  time  response  for  run 
number  11  in  which  R  **  10.  This  time  response  meets  the  required  specifications  for 
loi)  rate  Note,  however,  mat  the  gains  for  this  run  are  approximately  80%  higher,  on 

69 


the  average,  than  the  gains  for  the  best  run,  number  12,  in  Group  I.  in  order  to 
reduce  the  gains  so  that  only  a  small  control  effort  is  demanded,  there  is  more  work  yet 
to  be  dene. 

(3)  Group  3  Cost  Functions.  Making  the  penalty  on  the  Xj  error 
state  10  times  greater  than  the  penalty  on  the  and  Xj  error  states  in  the  Group  2 
cost  functions  appears  to  be  a  reasonable  mechanism  for  obtaining  an  adequate  time 
response.  In  an  effort  to  reduce  the  magnitude  of  the  control  effort,  the  ratios  of 

“^'^reased  to  100  for  the  Group  3  cost 
functions.  Table  8  and  Figure  3.13  present  the  data  for  this  group.  The  time  responses 
obtained  for  this  group  are  very  similar  to  those  obtained  for  Groups  1  and  2.  Notice 
in  Figure  3,13  that  the  overshoot  is  zero  for  all  cases  in  which  R  5.  In  addition, 
the  settling  time  is  less  than  one  second  when  R  £  20.  The  steady  state  gains  for  run 
number  14  average  only  42%  greater  than  for  run  number  12  in  Group  1.  Thus  the 
hypothesis  tested  in  the  Group  3  cost  functions  is  validated. 

(4)  Group  4  Cost  Functions.  The  cost  functions  in  Group  4  penalize  errors 
only  in  the  Xj  state  and  the  control  effort.  That  all  other  elements  of  H  and  Q  are  zero 
implies  that  no  penalty  is  assessed  against  deviations  in  the  X2  and  X3  states.  Table  9 
and  Figure  3.17  present  the  data  for  this  group.  Run  number  10  is  deemed  to  be  the 
most  acceptable  time  response  and  is  shown  in  Figure  3.19.  While  this  design  satisfies 
the  design  criteria,  note  that  the  steady  state  control  gains  average  93%  greater  than 
the  most  acceptable  run  in  Group  1.  This  is  the  greatest  increase  in  control  gains  that 
is  observed.  Also  notice  that  the  percent  overshoot  curve  in  Figure  3.17  inaeases 
upwards  of  9%.  This  is  much  greater  than  the  maximum  overshoot  of  4.65%  observed 
in  Groups  1,  2,  and  3.  For  these  two  reasons,  it  is  determined  that  the  cost  functions 
tested  in  Group  4  do  not  need  to  be  further  pursued. 

(5)  Summary.  Figures  3.21  and  3.22  summarize  the  information  contained 
in  Tables  6,  7,  8,  and  9.  It  is  interesting  to  note  in  Figure  3.21  that  the  first  three 
groups  of  cost  functions  yield  suprisingly  similar  curves  for  the  percent  overshoot  of 
the  roll  rate  system.  That  the  Group  4  cost  functions  produce  a  much  more  erratic 
curve  is  attributed  to  the  fact  that  no  weight  is  placed  on  the  Xj  or  Xj  states  in  this 
group.  The  roll  rate  settling  times  in  Figure  3.22  exhibit  similar  patterns  for  all  four 
groups  of  cost  functions.  Note  that  in  all  cases  there  appears  to  be  a  minimum  settling 
time  possible  when  2  ^  R  :S  10.  For  values  of  R  >  10,  the  large  emphasis  on  control 
effort  produces  small  steady  state  gains  which  in  turn  yield  a  slow  system. 
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Based  on  the  time  response  specifications  and  on  the  desire  to  design  a 
control  system  which  implements  minimal  steady  state  gains,  it  is  decided  that 
run  number  14  in  Group  3  is  the  best  solution  for  a  roll  rate  controller.  The  time 
response  for  this  set  of  parameters  appears  in  Figure  3.15. 

C.  AROD  ALTITUDE  RATE  CONTROLLER 
1.  The  Altitude  Rate  System 

•Because  the  primary  flight  mode  for  AROD  is  low  altitude  hover,  it  is 
important  that  there  be  a  reliable  control  system  to  maintain  the  vehicle's  vertical 
position  relative  to  the  earth's  surface.  The  throttle  on  AROD's  two  cycle  engine 
provides  the  mechanism  for  altitude  rate  control.  Table  10  defines  the  terras  which  are 
involved  in  the  altitude  rate  equations  of  motion. 


TABLE  10 

VARIABLE  DEFINITIONS  FOR  AROD  ALTITUDE  RATE  EQUATIONS 

OF  MOTION 


Variable 

Definition 

Value 

Units 

• 

h 

Vehicle  Altitude  Rate 

TBD 

feet/ second 

-Ch 

Engine  'Tlyrust  to.  RPM 
Dynamic  Coefficient 

0.0865 

ft/seconds^/rad 

Chapge  in 

Engme  Speed 

TBD 

RPM 

0.5 

seconds 

Ke 

Engine  Scale 

Factor 

837.8 

rad/ sec/rad 

• 

T^ottle  Servo 

Deflection  Angle 

^  30o 

radians 

». 

Throttle  Servo 

Deflection  Velocity 

S  50»/sec 

radians/second 

i 

TTirottle  Servo^  . 
Dampmg  Coefficient 

0.707 

unitless 

<0 

Throttle  Servo 

Natural  Frequency 

12.57 

radians/ second 

Control  Input 
to  Throttle  Servo 

TBD 

volts 
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By  commanding  a  desired  altitude  rate,  the  pilot  sets  in  motion  the  following 
sequence  of  events : 

1.  A  throttle  servo  control  signal,  u^,  is  generated  within  the  controller. 

2.  The  throttle  servo  position  is  adjusted. 

3.  The  actuator  position  commands  a  specific  engine  speed. 

4.  A  change  in  engine  RPM  causes  a  change  in  the  vehicle  altitude  rate. 

••  • 

The  rate  of  change,  h,  of  the  vehicle's  altitude  rate,  h,  is  proportional  to  the  change  in 
engine  RPM  as  shown  in  Equation  3.8. 

h  -  q  d,  (3.8) 

where  the  dynamic  constant,  is  experimentally  determined  in  wind  tunnel  tests. 
The  engine  is  modelled  as  a  first  order  lag  system  according  to  Equation  3.9. 

K  -  (- VV  ^  (Ke/te)  «t  (3-9) 

The  throttle  servo  is  modelled  as  a  second  order  plant  in  Equation  3.10. 

-  (o^6j  +  co^Uj  (3.10) 


The  signal  flow  graph  for  this  system  is  shovm  in  Figure  3.23.  The  following  state 
space  equations  are  used  to  design  the  controller  for  altitude  rate  ; 


(3.11) 
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0.0865 
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m  v 
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• 

0 

-2 

1675.5 
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0 

0 

1 

X  + 

0 

0 

0 

-157.91 

-17.77 
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mt  m 

(3.12) 


Figure  3.23  Signal  Flow  Diagram  for  Altitude  Rate  Control 
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(3.13) 


u 
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If  a  unit  step  is  commanded  for  the  altitude  rate,  then  the  command  input  vector 
becomes 


(3.14) 


2.  Altitude  Rate  Controller  Design 

The  system  given  in  Equations  3.12  and  3.14  is  entered  into  the  OPTCON 
program  and  a  controller  is  designed  according  to  a  procedure  similar  to  that  explained 
for  the  roll  rate  controller.  As  before,  a  sampling  rate  of  20  Hz  is  used  for  all  runs. 
Only  two  runs  are  hereafter  presented  because  the  lessons  learned  during  the  design  of 
the  roll  rate  controller  apply  equally  as  well  to  the  design  procedure  for  the  altitude 
rate  controller.  The  performance  specifications  for  this  control  system  are  designated 
to  be  as  follows  : 

1.  Zero  steady  state  error  for  a  step  input  is  required. 

2.  The  two  percent  settling  time,  t2o/^,  is  less  than  5  seconds. 

3.  No  overshoot  is  allowed. 

.  -The  first  run  is  made  using  a  baseline  cost  function.  The  results  obtained  for 
this  run  appear  in  Table  11.  Notice  that  an  incredibly  large  number  3f  stages  are 
required  in  order  for  to  be  achieved  using  this  cost  function.  If  the  gai)is  were  to  be 
implemented  dynamically  at  20  Hz,  more  than  38  seconds  would  be  required  before  the 
steady  state  gains  are  available.  This  clearly  is  not  desuable  since  the  settling  time 
must  be  less  than  five  seconds.  The  unit  step  time  response  using  stead}'  state  gains 
from  this  first  run  is  shown  in  Figure  3.24.  Even  after  20  seconds,  the  desired  altitude 
rate  is  not  yet  achieved.  The  cost  function  used  to  generate  this  solution  ii  deemed  to 
be  unsatisfactory  and  a  better  solution  is  sought. 

The  final  run  for  the  altitude  rate  controller  is  summarized  in  Table  12.  The 
cost  function  for  this  run  places  100  times  more  emphasis  on  errors  in  the  x^  state  than 
on  errors  in  the  X2  and  x^  states.  The  altitude  rate  time  response  shown  in  Figure  3.25 
exhibits  acceptable  performance  characteristics.  By  choosing  the  cost  function  wisely, 
it  becomes  possible  to  design  a  satisfactory  controller  for  this  system. 
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TABLE  1 1 

INITIAL  ALTITUDE  RATE  PARAMETERS 


Cost  Function 


“l 

1 

o 

o 

o 

“l 

0  0 

o' 

ft  • 

0 

1  0  0 

Q  - 

0 

1  0 

0 

R 

-  1 

0 

0  1  0 

0 

0  1 

0 

0 

0  0  1 

_0 

0  0 

1 

m 

Steady  State 
Gains 

Nunj^ber 

Percent 

Overshoot 

Settling 

Time 

fl 

^2 

^*3 

u 

Stages 

Required 

(sec) 

-0.0544 

-0.0461  -6.2776 

-0.1948 

763 

0.00 

>20.0 

TABLE  12 

FINAL  ALTITUDE  RATE  PARAMETERS 


Cost  Function 


100 

0 

0 

0 

m 

1 

0 

0 

««l 

0 

0 

1 

0 

0 

Q  - 

0 

.01 

0 

0 

0 

0 

1 

0 

0 

0 

.01 

0 

0 

0 

0 

1_ 

0 

0 

0 

.01  _ 

R  -  1 


Steady  State 
Gains 

^2  ^3 


Nuinber  Percent  Settling 
of  Overshoot  Tune 
Stages  ,  (sec) 

Required 


-0.3485  -0.0312  -4.5999  -0.1569 
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0.00 
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D.  AROD  PITCH  ANGLE  AND  YAW  ANGLE  CONTROLLER 
1.  The  Pitch  and  Yavt'  System 

Gyroscopic  coupling  between  the  pitch  and  yaw  dynamics  of  AROD  creates 
an  interesting  control  problem.  Refer  to  Table  13  for  an  explanation  of  the  terms 
involved  in  the  pitch  and  yaw  equations  which  follow. 


TABLE  13 

VARIABLE  DEFINITIONS  FOR  AROD  PITCH  AND  YAW 
EQUATIONS  OF  MOTION 


Variable 

Definition 

Value 

Units 

q 

Vehicle  Pitch  Rate 

TBD 

radians/second 

e 

Vehicle  Pitch  Angle 

TBD 

radians 

Pitch  to  Yaw 

-6.78 

seconds"^ 

q 

Gyroscopic  Coupling 

Elevator  Effectiveness 
Coefficient 

-14.51 

seconds*-^ 

• 

Eleyato.r  Servo . 
Deflection  Angle 

^  30-* 

radians 

«. 

Elevator  Seiyo  . 
Deflection  Veloaty 

S  50  • /sec 

radians/ second 

Control  Input 
to  Elevator  Servo 

TBD 

volts 

6 

r 

Vehicle  Yaw  Rate 

TBD 

radians/second 

V 

Vehicle  Yaw  Angle 

TBD 

radians 

Cr 

Yaw  to  Pitch 

Gyroscopic  Coupling 

6.75 

seconds’^ 
se  ;onds’^ 

K 

Rudder  Eftectiveness 

-16.68 

r 

Coefficient 

Rudder  Servo 

Deflection  Angle 

5  30o 

radians 

r 

• 

», 

Rudder  Servo 

Deflection  Velocity 

S  50®/sec 

radians/ second 

Control  Input 
to  Rudder  Servo 

TBD 

volts 

r 

; 

Elevator/ ladder  Servo 
Damping  Coefficient 

0.707 

unitless 

(1) 

EJcvatorJludder  Servo 
Natural  Frequency 

12.57 

radians/ second 

The  pitch  and  yaw  equations  of  motion  are  given  in  Equations  3.15  through  3.18. 


q  -  C^r  +• 


(3.15) 


e  -  J  q  dt 


(3.16) 


r  *•  C  q  +  N  6 
'-qM  r 


(3.17) 


V-frdt  (3.18) 

Note  that  crosscoupling  between  the  pitch  and  yaw  equations  enters  via  the  two 
gyroscopic  coupling  terms,  and  Cj,.  The  values  listed  in  Table  13  for  these  two 
coefficients  are  based  on  an  assumed  constant  propeller  velocity  of  72(X)  RPM. 

As  before,  the  elevator  and  rudder  control  vane  servos  arc  modelled  as  second 
order  systems  according  to  Equations  3. 19  and  3.20. 


5.  «  -2^(o5-  -  (0^6,  +  (o^u  (3.19) 

V  *9  w  '9 


dp  »  -2ija)6p  -  w^dp  +  (o^Up  (3.20) 


A  coupled  eighth  order  system  results  from  Equations  3.15  through  3.20.  The  signal 
flow  diagram  which  represents  this  MIMO  system  is  giver,  in  Figure  3.26. 
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Defining  the  eight  states  to  be  : 


X  -  q  5.  6,  V  r  «r  ^r]  ^ 

the  state  space  equation  for  the  pitch  and  yaw  coupled  system  is  defined  as 


(3.2!) 


X 


Ax  +  Ea 


(3.22) 


where 


0 

0 

0 

0 

0 

0 

0 

0 


''1 

0 

0 

0 

0 

-6,78 

0 

0 


0 

-14.01 

0 

-157.91 

0 
0 
0 
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0  0  0 

0  0  6.75 

i  G  0 

■  17,77  0  0 


0 

0 

0 
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0 

0 

0 

0 


-1 

0 

0 
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0 

0 

0 

0 

0 

-16.68 

0 

■  157.91 


0 
0 
0 
0 

0 

0 

1 

•17.77 


1 


(3.23) 


B 


0 

0 

0 

157.91 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

157.91 


(3.24) 


and  the  multi-input  control  vector  is 


-[i]- 


inuno 


{x-r} 


(3.25) 


2.  pitch  Angie  and  Yew  Angle  Controller  Design 
a.  Methodology 

Notice  in  Equation  3.25  that  the  control  input,  u,  is  a  (2  x  1)  vector.  Up 
to  this  point  in  the  AROD  control  design,  the  control  input  has  been  limited  to  a 
scalar  signal.  The  multi-input  control  that  results  from  gyroscopic  coupling  between 
pitch  and  yaw  requires  special  attention.  Consider  the  following  points : 

1.  The  optimal  feedback  gain  matrix,  F,  is  determined  from  the  solution  of  the 

discrete  matrix  Riccati  equation.  This  solution  requires  that  the  inverse  of  the 
4erm  P(K-l)  r  +  R}  in  Equation  2.28  be  determined.  .  - 

2.  The  cost  function  for  a  SI  SO  system  requires  that  the  control  weighting  factor, 
R,  be  a  scalar. 

3.  The  cost  function  for  an  n*  order  MI  MO  system  with  I  control  inputs  requires 
that  R  be  an  (C  x  J)  matrix. 

Thus,  for  a  SISO  system,  the  solution  for  F  is  greatly  simplified  because  the  term  in 
Equation  2.28  is  a  scalar  quantity.  For  a  MIMO  system,  however,  a  matrix  inversion 
routine  is  required  in  order  to  solve  for  the  optimal  gains.  Although  computationally 
possible,  it  is  decided  for  the  purpose  of  this  work  that  no  matrix  inversion  routine  is 
to  be  included  in  the  current  version  of  OPTCON.  This  implies  that  the  ability  of  the 
OPTCON  program  to  solve  for  optimal  feedback  gains  is  necessarily  limited  to  SISO 
systems.  This  limitation  is  reasonable  since  a  multitude  of  control  problems  can  be 
reduced  to  single  input  systems.  In  the  case  of  AROD,  however,  gyroscopic  coupling 
is  a  permanent  feature  of  the  pitch  and  yaw  dynamics.  Thus,  a  MIMO  system  is 
inevitable.  The  four  step  tactic  used  to  design  a  control  system  for  this  non-trivial 
problem  is  as  follows : 

1.  First  assume  that  the  gyroscopic  coupling  terms,  and  C^,,  are  zero  so  that  the 
coupled  eighth  order  system  reduces  to  two  independent  fourth  order  systems. 

2.  Use  OPTCON  to  solve  for  the  optimal  feedback  gains  for  the  two  independent 
systems. 

3.  Implement  the  steady  state  gains  so  obtained  for  the  fourth  order  uncoupled 
systems  into  a  simulation  model  for  the  eighth  order  coupled  system. 

4.  Experiment  with  various  combinations  and  modifications  to  the  (2  x  g) 
feedback  matrix,  Fj^^^q,  until  a  satisfactory  time  response  is  obtained  for  the 
pitch  angle  and  yaw  angle  of  the  coupled  system. 

Note  that  the  design  procedure  listed  above  do<*'  not  nt  the  most  direct  method 

to  design  a  MIMO  controller  using  optimal  contn.'  theory.  Rather,  this  method  is  an 

attempt  to  solve  a  complex  problem  using  a  tool  that  is  designed  to  solve  simpler 
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problems.  For  this  reason,  the  results  may  not  necessarily  be  expected  to  meet  the 
same  high  standards  required  of  the  two  previous  control  designs.  The  target 
performance  specifications  for  the  pitch  and  yaw  control  system  arc  stated  to  be  : 

1.  Zero  steady  state  error  for  a  step  input  is  required. 

2.  The  two  percent  settling  time,  t.2%>  ‘  seconds. 

3.  Less  than  10%  overshoot  is  allowed. 

(1)  Decoupling  the  System  Equations.  The  decoupling  procedure  results  in 
two  fourth  order  systems.  The  uncoupled  pitch  angle  state  space  equations  are  .  ^ 


X0  - 


(3.26) 


-1 

0 

0 

”  0  * 

• 

0 

0 

-14.51 

0 

0 

^  “ 

0 

0 

0 

1 

*0 

0 

0 

0 

-157.91 

-17.77 

157.91 

L  J 

(3.27) 


If  a  unit  step  is  commanded  for  the  pitch  angle,  then  the  pitch  command  input  vector 
becomes 


•’e 


(3.29) 
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The  uncoupled  yaw  angle  state  space  equations  are  : 


(3.30) 


“o 

-1 

0 

0 

0  1 

• .  • 

0 

0 

-16.68 

0 

0 

0 

0 

0 

1 

0 

0 

0 

-157.91 

-17.77^ 

J57.9lj 

(3.31) 


u. 


*’y) 


(3.32) 


If  a  unit  step  is  commanded  for  the  yaw  angle,  then  the  yaw  command  input  vector 
becomes 


(3.33) 


■The  similarity  between  the  dynamics  of  the  uncoupled  pitch  and  yaw  systems  is 
advantageous.  For  example,  it  is  found  that  the  elevator  and  rudder  control  gains, 
and  which  are  generated  by  OPTCON  have  identical  steady  state  values.  For  this 
reason,  only  one  set  of  steady  state  gains,  F,,,  needs  to  be  generated  for  each  cost 
function.  The  individual  elements  of  F^,  are  hereafter  referred  to  as  fj,  fj,  fj,  and  f^. 

(2)  Solving  for  the  Uncoupled  Controller.  The  solution  for  F,,  for  the  two 
fourth  order  systems  is  straightforward  and  follows  the  procedure  established  earlier  in 
this  chapter.  Table  14  summarizes  the  data  from  the  initial  run.  A  unit  step  time 
response  for  the  pitch  or  yaw  angle  controller  implementing  steady  state  gains  from 
Table  14  is  shown  in  Figure  3.27.  Table  15  contains  the  data  for  the  fmal  run.  The 
time  response  for  this  last  controller  appears  in  Figure  3.28.  The  steady  state  gains 
listed  in  Table  15  serve  as  the  foundation  upon  which  the  coupled  controller  is 
subsequently  designed. 


TABLE  14 

INITIAL  UNCOUPLED  PITCH  OR  YAW  PARAMETERS 


Cost  Function 


"l 

0 

0 

o" 

m 

1 

0 

0 

0* 

H  - 

0 

1 

0 

0 

Q  “ 

0 

1 

0 

0 

0 

0 

1 

0 

0 

0  . 

1 

0 

0 

0 

0 

1^ 

J) 
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0 

1^ 
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Stea^  State 
Gains 


Number 

or 

Stages 

Required 


Percent 

Overshoot 


Settling 

Time 

(sec) 


-0.1704  0.2417  -0.2490  -0.0858  96  0.00  4.15 


TABLE  15 

FINAL  UNCOUPLED  PITCH  OR  YAW  PARAMETERS 


Cost  Function 


H 


*100 

0 
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o“ 
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0 

o" 
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.01 
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5.0 


Steady  State 
Gams 
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I® 


equired 


loot  Time 
(sec) 


-0.3961  0.2808  -0.4158  -0.0259 
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Figure  3.27  Initial  Uncoupled  Pilch  or  Yaw  Angle  Time  Response 
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Figure  3.28  Final  Uncoupled  Pitch  or  Yaw  Angle  Time  Re$ponse 


TIME  (sec) 


(3)  The  Coupled  Feedback  Matrix.  In  order  to  implement  the  SI  SO 
feedback  gain  vectors,  and  F^,  into  the  coupled  MIMO  state  equations,  the  (2  x  3) 
feedback  matrix,  F  .  ,  is  formed  as  shown  in  Equation  3.34. 


(3.34) 


The  Djj  elements  of  represent  those  feedback  elements  which  are  not  specifically 
generated  by  OPTCON.  The  success  of  the  final  control  system  is  contingent  upon 
proper  selection  of  values  for  these  elements  of  F,^g.  For  the  purpose  of  the 
following  discussion,  the  reader  is  encouraged  to  refer  to  the  signal  flow  graph  shown 
in  Figure  3.26. 

(4)  Analysis.  There  are  numerous  ways  to  select  the  eight  unspecified 
feedback  gains  in  Equation  3.34.  The  seven  schemes  examined  during  the  course  of 
this  design  are  summarized  in  Table  16.  The  first  two  columns  in  Table  16  identify  the 
controller  structure  used  to  generate  the  elevator  and  rudder  control  signals,  u^  and  u^. 
The  third  column  refers  the  reader  to  the  appropriate  figure  containing  the  pitch  or 
yaw  angle  time  response  for  that  particular  set  of  parameters.  The  last  two  columns 
summarize  the  time  response  data  for  each  controller  design.  At  the  bottom  of 
Table  16  are  listed  the  exact  numerical  values  for  the  controller  gains. 
b.  Results 

(1)  Controller  Number  One.  The  feedback  matrix,  in  this  first 

design  assumes  that  the  four  states  of  the  yaw  equations  have  no  influence  on  the  pitch 
control  input,  uq.  Similarly,  the  four  states  of  the  pitch  equation  have  no  influence  on 
the  yaw  control  input,  Uy.  As  expected,  due  to  the  known  coupling  that  exists 
between  the  pitch  and  yaw  dynamics  of  the  vehicle,  the  time  response  in  Figure  3.29 
exhibits  unsatisfactory  performance. 

(2)  Controller  Number  Two.  For  this  design,  the  pitch  angle  state,  Xp 
influences  the  yaw  control  input,  Uy,  while  the  yaw  angle  state,  Xj,  contributes  to  the 
pitch  control  input,  uq.  From  the  time  response  in  Figure  3.30,  it  is  apparent  that  this 
design  is  not  satisfactory. 
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TABLE  16 

PITCH/YAW  CONTROLLER  DESIGN  SCHEMES 


Design 

Nuinber 

F 

tnimo 

Organization 

Time 

Response 

Figure 

Percent 

Overshoot 

S^yUng 

Time 

(sec) 

1 

‘3 

^■4 

0 

0 

0 

0 

3.29 

34.3 

18.9 

0 

0 

0 

0 

h 

f3 

^4 

2 

^"2 

^■3 

^■4 

h 

0 

0 

0 

3.30 

45.1 

N/A 

0 

0 

0 

h 

h 

^3 

U 

3 

h 

^*3 

C4 

h 

h 

0 

0 

3.31 

N/A 

N/A 

h 

0 

0 

h 

h 

h 

^4 

4 

h 

h 

^3 

^4 

fi 

h 

0 

0 

3.32 

2.15 

1.85 

h 

‘h 

0 

0 

h 

h 

h 

f4 

__  5 

h 

h 

^3 

Cj 

h 

h 

0 

0 

3.33 

0.00 

1.88 

h 

-f* 

0 

0 

h 

h 

f3 

^4 

6 

h 

^*2 

f3 

^4 

h 

^2 

0 

0 

3.34 

0.00 

N/A 

' 

h 

’h 

^3 

^4 

h 

h 

^3 

^4 

7 

h 

f3 

f4 

0 

h 

0 

0 

3.35 

22.5 

8.57 

0 

•^2 

0 

0 

^1 

^3 

f4 

fj  -  -0.3961 

1.2808 


Figure  3.29  Pitch  /  Yaw  Angle  Time  Response  for  Controller  Number  1 
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Figure  3.30  Pitch  /  Yaw  Angle  Time  Response  for  Controller  Number  2 
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Figure  3.32  Pitch  /  Yaw  Angle  Time  Response  for  Controller  Number  4 


95 


0.00  1.00  2.00  3.00  1.00  5.00 

REAL  TIME  (sec)  ' 


11  IX 


09‘I  02M  09’0  OIi’O  OO’O 

xaoi33rbyi  3ibis 


Figure  3.33  Pitch  /  Yaw  Angle  Time  Response  for  Controller  Number  5 
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Figure  3.34  Pitch  /  Yaw  Angle  Time  Response  for  Controller  Number  6 
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(3)  Controller  Number  Three.  Because  the  gyroscopic  coupling  between 
the  pitch  and  yaw  equations  is  directly  related  to  the  pitch  rate,  Xj,  and  the  yaw 
rate,  Xg,  it  is  decided  to  include  these  two  states  in  the  makeup  of  the  yaw  control  and 
pitch  control  respectively.  This  seems  like  a  logical  design  tactic  at  first  but  the 
resulting  time  response  in  Figure  3.31  proves  otherwise.  This  design  is  clearly  unstable. 

(4)  Controller  Number  Four.  Notice  in  the  coupled  system  signal  flow 

diagram  in  Figure  3.26  that  the  coupling  coeflicients,  and  C^,  are  nearly  equal  in 
magnitude  but  opposite  in  sign.  This  realization  causes  the  designer  to  hypothesize  that 
the  sign  of  Djj  in  should  be  reversed  from  the  value  previously  used  in 

controller  design  number  3.  As  shown  in  the  time  response  of  Figure  3.32,  this 
technique  yields  promising  results.  Although  a  steady  state  error  of  2.15%  exists,  there 
is  merit  in  pursuing  this  design  further. 

(5)  Controller  Number  Five.  By  finetuning  the  value  substituted  into  D22 
in  Fjjjjjjjg,  the  time  response  for  pitch  angle  or  yaw  angle  is  made  to  satisfy  the  desired 
performance  criteria.  The  time  response  in  Figure  3.33  exhibits  no  steady  state  error, 
zero  overshoot,  and  a  settling  time,  tjo/^,  of  slightly  less  than  two  seconds.  This 
controller  design,  then,  is  completely  satisfactory. 

(6)  Controller  Number  .Six.  For  this  design,  a// states' are  allowed  to 
influence  both  u^  and  u^.  The  time  response  so  obtained  is  shown  in  Figure  3.34. 
Even'  though  the  steady  state  angle  is  only  75%  of  the  commanded  value,  this 
controller  design  appears  to  be  potentially  useful.  By  tuning  the  gains  iteratively,  it  is 
believed  that  zero  steady  state  error  is  achievable  with  this  design. 

(7)  Controller  Number  Seven.  This  final  design  is  a  modification  to 
controller  number  3.  In  this  case,  only  the  pitch  rate  and  yaw  rate  contribute  to  the 
crosscoupled  control  signals.  This  design  effort  results  in  unsatisfactory  performance 
as  shown  in  Figure  3.35. 

c.  The  Final  Design 

Controller  number  5  is  selected  as  the  best  design  for  a  pitch/yaw  angle 
controller.  The  time  response  for  this  design  appears  in  Figure  3.33.  Note  that  this 
design  is  based  on  feedback  gains  generated  by  OPTCON  but  that  a  modification  to  f^ 
is  required  in  order  to  obtain  the  final  design.  Thus,  the  controller  is  not  optimal,  by 
formal  definition,  even  though  optimal  control  theory  provides  the  foundation  for  its 
development. 
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IV.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

The  lessons  learned  during  the  course  of  this  work  are  as  follows  : 

1.  The  cost  function  weighting  parameters,  H,  Q,  and  R,  play  vital  roles  in 
determining  the  magnitude  of  the  steady  state  optimal  feedback  gain  matrix, 
Fjj.  These  control  gains,  in  turn,  significantly  affect  the  time  response- of  the 
controlled  system. 

2.  There  is  no  magic  formula  to  determine  proper  values  for  the  weighting  factors. 
A  reasonable  starting  point  is  to  use  the  baseline  cost  function  in  which  all 
diagonal  elements  of  H,  Q,  and  R  are  assigned  the  value  of  unity  and  all  off- 
diagonal  elements  are  zero. 

3.  The  process  of  trial  and  error  is  prerequisite  to  the  successful  design  of  an 
optimal  control  system.  Only  through  an  iterative  procedure  does  the  engineer 
establish  the  true  nature  of  the  problem  at  hand. 

4.  There  are  obvious  trends  to  be  aware  of.  These  include  : 

a.  The  sampling  frequency,  f,,  must  be  fast  enough  to  avoid  aliasing  effects. 
As  predicted,  the  use  of  a  sampling  frequency  that  is  five  to  ten  times  faster 
than  the  Nyquist  frequency  seems  to  be  adequate. 

b.  As  the  selected  sampling  frequency  increases,  the  optimal  gains  generated 
also  tend  to  increase  in  magnitude. 

c.  The  control  weighting  factor,  R,  for  a  SISO  system  can  be  used  as  a 
parameter  to  systematically  alter  the  time  response  of  the  system.  As  the 
relative  weight  on  the  control  effort  increases,  the  steady  state  gains  tend  to 
decrease  in  magnitude.  This  generally  produces  a  slow  system  that  exhibits 
little  or  no  overshoot.  On  the  other  hand,  if  the  value  of  R  is  decreased, 
the  steady  state  gains  can  become  so  large  that  a  very  fast  and  highly 
oscillatory  system  results. 

5.  The  controller  design  for  a  MIMO  system  is  significantly  more  involved  than 
the  design  for  a  SISO  system.  If  the  engineer  can  logicaUy  and  accurately 
decouple  the  MIMO  system  into  multiple  SISO  systems,  then  the  design  effort 
becomes  much  easier.  As  shown  in  the  pitch  and  yaw  controller  for  AROD, 
this  method  is  feasible. 

B.  RECOMMENDATIONS  FOR  FURTHER  WORK 

The  following  areas  present  valid  opportunities  for  useful  expansion  of  this 
work  : 

1.  A  parameter  identification  procedure  which  aids  the  design  engineer  in 
determining  or  estimating  the  A,  B,  <P,  and  F  plant  matrices  is  needed.  The  use 
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of  a  Fast  Fourier  Transform  (FFT)  algorithm  is  one  possible  solution  to  this 
requirement.  Such  a  program  could  be  used  in  conjunction  with,  but  not 
necessarily  integrated  into,  the  existing  OPTCON  package. 

2.  The  OPTCON  program  is  limited  in  that  it  does  not  generate  optimal  feedback 
gains  for  a  .MI MO  system.  A  matrix  inversion  routine  is  needed  so  that  the 
discrete  matrix  Riccati  solution  can  be  determined  for  any  (n  x  £)  B  or  <D 
matrix. 

3.  The  theory  of  optimal  control  assumes  availability  of  all  states  for  feedback. 
The  design  process  must  account  for  the  fact  that  all  states  are  not  always 
measured.  In  the  case  of  AROD,  the  servo  position  and  rate  states  are  not 
available  for  feedback.  This  means  that  an  observer  must  be  designed  iH  order 
to  provide  the  missing  state  information. 

4.  The  three  control  systems  which  are  herein  designed  must  be  evaluated  on  the 
actual  vehicle.  Although  computer  simulations  provide  a  wealth  of  insight,  the 
proof  of  a  good  design  rests  in  the  ability  of  the  system  to  function  in  the 
outside  world. 
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APPENDIX  A 
THE  OPTCON  PROGRAM 


1.  OVERVIEW 

The  purpose  of  this  appendix  is  to  describe  in  detail  the  OPTCON  computer 
program  which  was  developed  in  support  of  this  thesis.  OPTCON  derives  its  name 
from  OPTimal  CONtrol.  A  previous  edition  of  OPTCON  by  Professor  H.A.  Titus  of 
the  Naval  Postgraduate  School  provided  the  starting  point  for  the  work  that  follows. 
The  original  OPTCON  program  allowed  the  user  to  input  a  state  space  system  either  in 
the  continuous  time  domain  or  in  the  discrete  time  domain.  Using  matrix  calculations 
to  solve  Equations  2.28,  2.29,  and  2.30,  this  first  version  of  OPTCON  generated  a  table 
of  feedback  gains  and  immediately  terminated  execution.  The  motivation  for 
improving  the  original  OPTCON  is  fourfold. 

1.  The  design  process  is  an  iterative  technique.  The  OPTCON  program  needs  to 
be  flexible  enough  to  allow  minor  changes  to  be  made  to  specific  parameters 
without  the  requirement  to  re*initialize  all  of  the  cost  function  and  system 
values. 

2.  The  gain  trajectory  table  is  not  a  convenient  means  by  which  to  analyze  the 
solution  to  an  optimal  control  problem.  A  graph  of  the  feedback  gains  verses 
time  provides  better  insight. 

3.  A  time  response  of  the  state  space  is  needed  in  order  to  allow  the  designer  to 
quickly  evaluate  the  performance  of  the  system. 

4.  The  program  should  be  user  friendly.  The  original  OPTCON  demanded  that 
the  user  flawlessly  enter  the  correct  response  to  every  question  on  the  first 
attempt.  Woe  be  it  to  the  user  who  accidentally  types  a  letter  in  response  to  a 
question  that  requires  a  numerical  answer.  The  program  aborts  and  any  effort 
that  was  spent  in  entering  information  is  wasted.  The  frustration  factor  for 
such  an  unfriendly  program  is  likely  to  leave  the  program  sitting  on  the  shelf 
with  nobody  to  use  it. 

With  these  points  in  mind,  the  OPTCON  program  is  rewritten  to  provide  an 
interactive,  menu  driven,  user-friendly,  optimal  control  design  tool  that  capitalizes  on 
the  graphical  capabilities  of  modern  microcomputers.  The  program  is  written  in 
MICROSOFT  Fortran  and  is  listed  in  Appendices  B,  C,  and  D.  Appendix  B  contains 
the  driver  program,  MAIN.  Appendix  C  contains  the  majority  of  the  subroutines 
which  are  called  by  MAIN.  Appendix  D  contains  the  plotting  subroutine,  GRAPH, 
which  makes  use  of  the  PLOT88  graphics  software  package. 
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In  order  to  use  OPTCON  to  its  full  potential,-  the  user  needs  access  to  the 
following  : 

1.  A  microcomputer  capable  of  executing  MICROSOFT  Fortran  based  programs. 
During  the  development  of  OPTCON,  an  IBM  AT  computer  configured  with 
640  kbyte  RAM  and  Intel's  80287  math  co-processor  was  used. 

2.  Fortran,  PLOT88,  and  Math  libraries. 

3.  A  monochrome  or  color  graphics  monitor. 

4.  An  Epson  or  LaserJet  printer. 

Figure  A.l  is  provided  to  give  the  user  a  broad  overview  of  the  basic  program  -flow  of 
OPTCON.  The  blocks  outlined  by  solid  lines  represent  program  segments  that  must  be 
performed  during  the  initial  execution  of  OPTCON.  The  blocks  outlined  by  dashed 
lines  represent  program  segments  that  are  optional.  The  numbers  that  appear  to  the 
leR  of  each  block  are  referred  to  during  in  Section  2.d  of  this  appendix. 

The  remainder  of  this  appendix  illustrates  the  solution  to  a  simple  example 
problem  using  the  OPICON  program.  The  intent  here  is  not  to  focus  on  the  specific 
example  problem  or  on  its  solution  but,  rather,  to  focus  on  the  capabilities  and  use  of 
OPTCON. 

2.  AN  EXAMPLE  PROBLEM 
a.  The  Second  Order  Integrator 

Consider  the  continuous  time  system  shown  in  Figure  A.2.  The  state  space 
equations  for  this  second  order  plant  are  derived  by  defining  the  output  of  each 
integrator  to  be  a  state.  Using  Equations  2.7  through  2.10  as  a  basis,  the  state 
equations  become  : 


^'>■[0  0]^ 

-[0  3 


X(t)  u(t) 

(A.1) 

X(t) 

(A.2) 

■1 

+ 

1 

(A.3) 
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Figure  A.l  OPTCON  Program  Flow  Diagram 
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_  Figure  A.2  Second  Order  Integrator  Signal  Flow  Diagram 
If  the  system  is  sampled  every  T  seconds.  Equation  2.20  yields : 


and  Equations  2.21  and  2.22  yield 


(A.4) 


(A.5) 


(A.6) 
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The  preceeding  calculations  are  done  simply  to  allow  verification  of  the  PHIDEL 
subroutine  in  Appendix  C.  This  subroutine  converts  an  A  and  B  continuous  system  to 
a  O  and  F  discrete  system.  For  instance,  assume  that  the  system  is  sampled  at 
fj  =  100  Hertz.  This  means  that  T  =  0.01.  Equations  A.5  and  A.6  become 


for  the  discrete  time  representation  of  the  second  order  integrator. 

Before  proceeding  with  the  OPTCON  program,  the  user  is  urged  to  verify  that 
the  system  is  controllable  and  observable, 
b.  Controllability  and  Reachability 

According  to  Astrom,  a  system  is  controllable  [Ref.  5:  p.  104]  only  if  "it  is 
possible  to  find  a  control  sequence  such  that  the  origin  can  be  reached  from  any  initial 
state  in  finite  time."  Thus,  controllability  is  a  necessary  condition  for  the  regulator 
problenx  A  similar  property  called  reachability  is  required  for  the  tracking  problem.  A 
system  is  reachable  [Ref.  5:  p.  104]  only  if  "it  is  possible  to  find  a  control  sequence 
such  that  an  arbitrary  state  can  be  reached  from  any  initial  state  in  finite  time." 

— For  continuous  time  systems,  the  properties  of  controllability  and  reachability 
coincide.  That  is,  cither  a  continuous  time  system  is  both  controllable  and  reachable 
or  it  is  both  uncontrollable  and  unreachable.  For  a  discrete  time  system,  however, 
controllability  does  not  guarantee  reachability.  Reachability  of  a  discrete  time  system 
does  guarantee  controllability.  The  reachability  of  a  discrete  system  is  important 
because  the  engineer  should  not  spend  a  lot  of  time  designing  a  controller  that  is 
physically  impossible  to  implement. 

A  simple  test  is  performed  to  check  the  reachability  of  a  discrete  system.  The 
(n  X  nf)  reachability  matrix,  W^.,  for  an  n*  order  discrete  time  system  with  t  control 
inputs  is  defined  as  follows  : 

-  [jr  or  o^r  ...  (a.9) 
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If  th?  reachability  matrix  is  of  rank  n,  then  the  system  is  reachable.  In  the  case  of  the 
example  problem 


(A.  10) 


Taking  the  determinant  of  W^.  and  setting  the  result  equal  to  zero,  the  necessary 
condition  for  reachability  is  found  to  be  that  T  *  0.  Since  an  infinite  sampling 
frequency  is  impossible  to  achieve,  the  system  is  reachable  and  it  is'  reasonable  to 
continue  with  the  problem, 
c.  Observability 

In  order  to  take  full  advantage  of  the  optimal  gain  schedule,  F(k),  it  is 
necessary  that  all  of  the  states  be  observable.  According  to  VanLandingham 
[Ref.  4:  p.  308],  a  discrete  time  system  is  completely  observable  if  it  is  possible  to 
determine  the  initial  state,  x(0),  of  the  system  based  on  knowledge  of  the  control  input, 
u(k),  and  the  output,  y(k),  over  a  finite  number  of  time  intervals.  The  test  for 
observability  closely  follows  the  test  for  reachability.  First,  define  the  (mn  x  n) 
observability  matrix,  W^,  as 


Wo 


c 

CO 

co^ 


(A.  11) 


If  the  rank  of  is  n,  then  the  system  is  observable.  This  implies  that  all  of  the 
states  of  the  system  are  available  for  state  feedback.  If  the  system  is  not  completely 
observable,  then  one  or  more  of  its  states  is  not  measureable.  Either  the  system  must 
operate  without  the  unobserved  states  in  the  feedback  path,  or  an  observer  must  be 
designed  to  estimate  the  unobserved  states.  In  the  example  problem,  the  observability 
matrix  is 


W. 


I  0 
0  1 
1  T 

0  1 
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(A.  12) 


107 


This  observability  matrix  is  of  rank  2  and  the  system  is  completely  observable.  Notice 
that  if  the  output  distribution  matrix  is 


C 


(A.13) 


so  that  only  Xj  is  observed,  the  observability  matrix  becomes 


which  is  of  rank  2  provided  that  T  * 

C  -  [o 


and  the  system  is  not  observable  regardless  of  the  sampling  frequency.  A  state 
observer  is  needed  to  give  an  estimate  of  the  Xj  state. 

d.  Solution  Using  OPTCON 

_This  section  is  an  introductory  guide  to  OPTCON.  The  second  order 
integrator  problem  is  used  to  acquaint  the  user  with  the  commands,  features,  and 
limitations  of  the  program.  The  messages  presented  in  this  section  are  referred  to  as 
"screens"  and  are  surrounded  by  numbered  boxes.  Neither  the  boxes  nor  the  numbers 
by  which  they  are  referenced  are  actual  features  of  the  OPTCON  program.  They  are 
simply  used  as  devices  to  make  the  following  discussion  more  understandable. 
Messages  which  are  generated  by  OPTCON  appear  in  standard  print.  Any  responses 
which  represent  keyboard  entry  by  the  user  are  shown  in  italic  print.  If  the  re^iponse  is 
to  be  y  for  "yes"  or  n  for  "no",  then  either  uppercase  or  lowercase  letters  are  acceptable. 
If  the  response  is  to  be  an  integer  entry,  as  in  the  menu  selections,  the  subprogram 
COMPARE  is  called  to  verify  that  the  user  has  entered  a  valid  integer.  If  the  response 
is  out  of  range  of  the  acceptable  values,  or  if  the  response  is  not  an  integer,  then  the 
program  repeats  the  message  until  a  valid  response  is  entered  by  the  user. 


0.  However,  if  only  Xj  is  observed,  then 


0 

;] 


(A.15) 

(A.  16) 
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/.  Starting  the  Program 

The  user  enters  OPTCON  by  typing  optcon  on  the  command  line.  The 
following  heading  appears 


Screen  1 

OPTCON  ainlaiaM  Vw  folXoMing  ocm\  funotloni 

J  alflN  (  XMN)  «  H  •  X(N)  *  Sumt  X‘(K)  «  0  »  X(k)  *  UMK)  «  R  «  U(k))i 

Th«  output  of  tho  prograa  is  tho  foodbooK  gain  MitriP*  P  tronopooa*  (P'}> 
Hhieht  Mhon  ■ultipllod  by  tho  Stato  Vaetor  IX)» 
yialda  a  scalar  eontrel.(U). 

Tha  folloMlng  raeursiva  aquations  Mara  dsrivad  using  dynaalo  pregraoMingy 
starting  at  tha  taralnal  tiao  (N)  and  Marking  baokMordst 

(1)  P'(R)  ■  -(OEL*«P(k-l)«PHZ)/(OeL*«P(k-I)<iOEL  ♦  R)  PM0)a0 

(2)  PSl(k)  >  PHZ  ♦  OELaPMk)  PSZ(0)a0 

(S)  P(k)  >  PSZ‘(kJ«P(k-l)«PSZ(k)  *  a  *  FMk)«RaP(k)  P(0)>H 


2.  Entering  Initial  htformation 

The  first  entry  required  is  a  problem  name.  This  name  is  used  to  identify 
the  output  file  called  OPTFILE  which  contains  all  matrices,  gain  trajectories,  and  time 
response  trajectories  for  each  run  that  the  user  requests  during  the  problem  session. 


Screen  2 

First  antar  tho  preblaa  idantifioation  (  NOT  to  aicoasd  20  ohoroetars  ). 
PROBLEN  10........ second  order  example 


Next,  the  user  selects  the  type  of  printer  that  is  connected  to  the  operating 
system.  If  graph  hardcopies  are  not  to  be  requested  during  the  course  of  the  problem 
session,  then  the  response  to  this  question  has  no  impact  on  the  operation  of 
OPTCON.  If  graph  hardcopies  are  to  be  requested,  however,  then  the  answer  to  this 
question  sets  a  flag  that  allows  data  to  be  properly  formatted  for  the  printer  that  is 
being  used.  Unpredictable  results  are  expected  if  the  user  attempts  to  get  graph 
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hardcopies  from  a  printer  that  is  not  selected.  In  this  example,  the  LaserJet  printer  is 
to  be  selected. 


Now  the  user  enters  the  order,  n,  of  the  system.  The  maximum  order 
which  OPTCON  can  accept  is  eight  due  to  the  limitation  of  64  kbytes  per  segment  in 
the  IBM  AT  microcomputer.  The  practice  problem  requires  that  a  2  be  entered  here. 


3.  Entering  the  Cost  Function 

Next,  the  number  of  discrete  time  stages,  N,  is  entered.  This  number  is 
limited  TO  1000  due  to  the  maximum  dimension  size  of  the  arrays  in  OPTCON.  The 
user  should  be  aware  of  the  relationship  : 

tf  -  NAt  (A.  18) 

where 

tf  -  final  time  of  the  process 

N  "  number  of  stages  over  which  the  T  in  Equation  2.23 
IS  to  be  performed. 

At  *  sampling  interval 

In  the  example,  let  At  ■  0.01  seconds  and  tj.  10  seconds.  This  requires  N  ■  1000. 


I 

*►« 


At  this  point,  the  weighting  elements  of  the  cost  function  are  to  be  entered. 
Assuming  that  the  user  wants  to  initially  create  a  baseline  solution  for  the  problem,  a 
reasonable  starting  point  is  to  let  all  diagonal  weighting  factors  assume  a  value  of 
unity.  The  routine  to  enter  the  cost  function  begins  with  Screen  6. 


Screen  6 

0o«s  oest  function  (J)  include  the  Stete  TRAJECTONY  over  ell  etegee  ? 

(  Anener  l>2>er  3  ) 

1)  YES... Set  a  equel  to  the  IDEMTXTY  Hetrix  . 

2)  YES.. .Eeoh  dieflonel  eloeent  of  9  will  be  entered  eeperetely. 

S)  NO _ Set  a  eciuBl  to  the  ZERO  Hetrix  . 


ANSNER. 


Selecting  option  1  results  in  th<s  Q  matrix  being  echoed  in  Screen  7.  The 
program  then  advances  directly  to  Screen  11. 


Screen  7 

The  etetee  are  weighted  equally  for  the  TRAJECTORY  over  all  etegee. 

The  a  Hetrix 

1.0000  .0000 
.0000  1.0000 


Selecting  option  2  in  response  to  Screen  6  allows  the  user  to  enter  values 
for  the  diagonal  elements  of  the  Q  matrix.  All  ofT-diagonal  elements  are  automatically 
set  equal  to  zero.  For  the  sample  problem,  assume  that  the  user  wants  ■  2.4  and 
q22  “  5.  After  entering  the  value  for  q^,  the  user  is  prompted  to  enter  the  value  for 
q22-  Screen  8  results. 


,^L 
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Screen  8 


entap  'tha  alaaanlt  of  tha  t)  aatrix. 

(Sfafa  Maighfing  aatrlx  for  TRAJICTQftY  evar  all  afagaa) 

e(l>ll  >  T  2.4 
.9(2,2)  >  ?  5 


Afler  the  user  enters  all  diagonal  elements,  the  matrix  is  echoed  in  Screen  9. 
OPTCON  then  advances  to  Screen  11. 

Screen  9 

Tha  9  Hafpix 

2.4000  .0000 

.0000  5.0000 


Selecting  option  3  in  response  to  Screen  6  sets  all  elements  of  the  Q  matrix 
equal  to  zero.  Screens  10  and  1 1  follow. 

Screen  10 

Tha  sfa4a  TfUUECTORY  is  not  Inoludad  in  your  ooat  fuwtlon. 

Tha  0  Matrix 

.0000  .0000 

.0000  .0000 


Screen  1 1  involves  a  loop  which  allows  the  user  to  change  any  or  ail  of  the 
elements  in  the  matrix  that  is  currently  being  processed.  This  loop  is  subsequently 
referred  to  in  this  discussion  as  "the  modify  routine." 


Screen  1 1 

Do  you  Mont  to  ohangs  any  alaaani  of  too  rntrlxt 

.1)  Yes...a  aZNOLB  aloMnt. 

2)  YIS...«ha  INTZRI  HatolK. 

3)  NO 

ANSMa . 

Option  1  produces  Screen  12  which  allows  the  user  to  change  a  single 
element  by  identifying  the  row  and  column  of  the  element  to  be  changed.  The  row  and 
column  entries  must  be  integers  separated  by  a  comma.  Assume  that  the  user  wants  to 
change  q22  so  that  it  equals  3. 


«  Screen  12 

alaMant  of  fho  rix  do  you  Mont  fo  Ction0o  ? 

Zf  Z  ia  tho  RON  and  J  is  ttw  COLUMN . .  Z>J  2,2 


The  user  is  then  prompted  to  enter  the  new  value  for  the  element  that  is  to 
be  changed.  Screen  13  applies. 


Screen  13 


ef2»Z)  >73 


If  the  user  entered  this  situation  directly  from  Screen  7  then  the 
result  is  Screen  14. 


Screen  14 

Tha  9  Matrix 

1.0000  .0000 

.0000  3.0000 

Any  othar  ohangaa?  (AnaMar  y  or  n) 

1*  . 

- 

If  the  answer  to  Screen  14  is  j/,  then  OPTCON  returns  to  Screen  12  and 
allows  changes  to  be  made  to  individual  elements  of  the  matrix.  Once  the  user  is 
satisfied  that  that  the  Q  matrix  is  correct,  a  n  is  entered  in  response  to  Screen  14.  At 
this  point,  OPTCON  is  ready  to  accept  information  relating  to  the  H  matrix. 
Screen  15  is  next. 


Since  this  is  the  desired  H  matrix,  &  n  is  entered  and  the  program  advances 
to  the  section  in  which  the  user  is  asked  to  enter  the  value  for  R.  Assume  the  desired 
value  is  to  be  15.7.  Screens  17  and  18  result. 


The  program  echoes  the  value  entered  and  asks  if  there  are  any  changes.  If 
there  are  changes  to  be  made,  ay  response  returns  the  user  to  Screen  17.  An  response 
in  Screen  18  indicates  that  the  cost  function,  J,  is  now  defined  completely  and  Block  1 
of  Figure  A.l  is  complete.  The  program  advances  to  Block  2  of  Figure  A.l 

The  user  must  now  indicate  if  the  problem  to  be  solved  is  in  the  continuous 
time  d.0Qiain  or  in  the  discrete  time  domain.  Screen  19  applies. 


The  sample  problem  is  of  the  continuous  type  and  a  0  is  the  appropriate 
response  to  Screen  19.  Screen  20  follows. 
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If  a  response  is  entered  in  Screen  20,  then  the  program  advances  to 
Screen  22.  Otherwise,  the  message  in  Screen  21  appears. 


Screen  21 


You  will  anUr  iha  PHX  and  DEL  aairioM. 
. It  this  oorraot  7 


The  program  toggles  between  Screens  20  and  21  until  the  user  enters  y  to 
one  of  these  two  options.  Assuming  that  the  continuous  system  is  selected,  the  next 
section  of  the  program  allows  the  user  to  enter  the  A  and  B  system  matrices  and  the 
sampling  interval,  T. 

4,  Entering  the  Continuous  Time  System  Parameters 

The  elements  of  the  A  matrix  are  sequentially  entered  as  shown 
in  Screen  22. 


Screen  22 

Entsr  th«  wImm 

Hits  of  fho  plant  Matrix— A. 

A(l>ll  «  7 

0 

• 

Ad, 2)  ■  7 

1 

A(2,l)  ■  7 

0 

A(2,2I  a  7 

0 

Screen  23  echoes  the  A  matrix  and  affords  the  user  an  opportunity  to  make 
any  changes.  The  modify  routine  is  entered  unless  the  user  responds  to  Screen  23  with 
a  3.  In  the  case  shown,  all  entries  are  correct  and  a  J  is  appropriate. 
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Screen  23 

Tha  A  Matrix  (Plant  Matrix  1 

.0000  1.0000 

.0000  .0000 

Do  you  Ma.'it  to  ohanea  any  alanant  of  tha  Matrix? 

11  YES... a  SZNQLE  alanant. 

2 1  YES. . . tha  ENTIRE  Matrix. 

SI  NO 

ANSMBR . J 

The  elements  of  the  B  matrix  are  sequentially  entered  as  shown 
in  Screen  24. 


Screen  24 

Entar  tha  alaai 

anta  of  tha  diatribution  natrix-*-B. 

B(ia)  •to 

B(l»l)  •  t  i 

Screen  25  echoes  the  B  matrix  and  once  again  allows  the  user  to  enter  the 
modify  routine  if  necessary. 


Screen  25 


Th«  B  tlitriN  ( Distrlbu'tian  Hitrixl 

.OOCO 

1.0000 

Do  you  Marti  to  ohanea  any  alaaani  of  iha  Matrix? 
II  Ye8...a  SZNBLI  alaMant. 

21  Yes...tha  enrzu  Matrix. 

SI  NO 


ANSNER. 


Since  no  changes  are  needed,  the  program  now  prompts  the  user  to  enter 
the  sampling  interval,  T.  The  correct  answer  for  the  sample  problem  is  entered 
in  Screen  26. 
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As  usual,  the  response  is  echoed  and  the  user  is  allowed  to  make  changes 
until  the  correct  value  is  entered.  Screen  27  applies. 

A  . 

Screen  27 

Th«  SAMPLE  INTERVAL  DT  >  .0100 

Any  Chang**  to  th*  SAMPLE  INTERVAL  7  (An*M*r  y  or  nil  n 

5,  The  Optimal  Feedback  Gains  Calculated 

The  program  now  has  all  of  the  information  necessary  to  calculate  the 
optimal  gain  schedule.  The  first  step  that  OPTCON  must  perform  is  to  convert 
the  A  and  B  matrices  to  the  corresponding  0  and  F  matrices.  The  subroutine 
PHIDEL  in  Appendix  B  performs  this  conversion.  The  resulting  O  and  F  matrices  are 
not  displayed  on  the  monitor.  These  two  matrices  are,  however,  recorded  in  the 
OPTFILE  listing  for  the  user's  convenience.  If  a  discrete  time  system  is  initially 
selected  in  Screens  19  and  20,  then  the  PHIDEL  subroutine  is  bypassed.  In  either 
case,  the  gain  schedule  is  now  calculated  using  Equations  2.28,  2.29,  and  2.30.  This 
completes  Block  3  of  Figure  A.l.  As  block  4  of  Figure  A.l  is  entered,  the  user  may 
choose  to  view  the  gain  schedule  in  tabular  fomi.  Screen  28  applies. 

Screen  28 

Do  you  Mant  to  a**  th*  gain  aehadula  tabl*  on  tha  aeraan  7 
(AnaMar  y  or  n)  y 

Since  the  user  wishes  to  view  the  gain  schedule  table  on  the  monitor,  a  is 
entered  in  Screen  28.  The  user  should  remember  two  points  before  choosing  to  list  the 
gain  schedule  on  the  screen  : 

1.  The  gain  schedule  is  automatically  entered  into  OPTFILE  regardless  of  the 
user's  response  in  Screen  28.  If  the  user  wants  to  record  the  exact  values  of  the 


gains,  this  output  file  may  be  examined  later  using  the  BROWSE,  COPY, 
EDIT,  PRINT,  or  TYPE  commands  in  DOS. 

2.  A  total  of  N  lines  of  data  are  listed  on  the  monitor  when  tabular  output  is 
requested  in  Screen  28.  If  N  is  on  the  order  of  several  hundred  or  more,  the 
design  procedure  is  likely  to  lose  momentum  due  to  the  lengthy  delay  involved 
in  sending  such  a  large  array  to  the  monitor. 

In  order  to  illustrate  the  form  of  the  data  generated.  Screen  29  lists  a 
portion  of  the  gain  schedule  table.  Only  the  first  ten  time  intervals  are  listed  here  for 
brevity.  The  actual  sample  problem  lists  a  table  with  1000  rows. 


Screen  29 


NEO 

TIME 

STEP 

REAL 

TIME 

INDEX 

FMl) 

FM2) 

1 

1000 

.0000 

-.0100 

2 

999 

-.0002 

-.0200 

3 

998 

-.0004 

-.0300 

4 

997 

-.0008 

-.0400 

5 

994 

-.0012 

-.0500 

4 

99S 

-.0018 

-.0400 

7 

994 

-.0024 

-.0700 

8 

993 

-.0032 

-.0800 

9 

992 

-.0040 

-.0900 

10 

991 

-.0050 

-.1000 

Block  4  of  Figure  A.l  is  now  complete  and  Block  5  is  initiated.  The  next 
option  available  to  the  user  is  to  have  OPTCON  generate  graphs  of  the  optimal  gain 
trajectories.  If  graphs  are  not  desired,  the  user  may  answer  n  in  response  to  Screen  30 
and  the  progr^  advances  to  Screen  32.  If  plots  of  the  gain  trajectories  are  desired, 
then  a  y  response  is  required  in  Screen  30. 


Screen  30 


Do  you  MWit  to  tha  gwina  plettad  7 
(Anawar  y  or  nl  y 


At  this  point,  the  program  calls  subroutine  GRAPH.  An  internal  flag  is  set 
so  that  the  gain  trajectory  plots  are  sent  to  the  monitor.  A  separate  plot  is  generated 
for  each  gain  trajectory.  Thus,  for  an  n*  order  system,  there  are  n  separate  gain  plots 
produced.  As  each  graph  is  generated  on  the  screen,  a  pause  is  inserted  so  that  the 
user  may  conveniently  examine  each  one.  Striking  any  key  removes  the  current  graph 
from  the  monitor  and  Screen  3 1  follows. 

Screen  31 

Do  you  Nant  ■  hardoopy  of  fhia  plot  ?  C  AnaMor  y  or  n  )  y 

If  a  /I  is  entered  in  Screen  31,  then  the  program  begins  to  generate  the  next 
gain  trajectory  plot  for  monitor  output.  By  answering  y  in  Screen  31,  the  user  will 
automatically  be  provided  with  a  hardcopy  of  the  gain  trajectory.  A  single  hardcopy 
graph  requires  approximately  120  seconds  on  the  Epson  printer  and  approximately  90 
seconds  on  the  Laserjet  printer.  Because  of  the  superior  quality  of  the  graphs  available 
from  the  later,  all  graphs  contained  in  this  thesis  are  generated  on  the  Laserjet  printer. 
As  soon  as  the  hardcopy  is  complete,  OPTCON  begins  to  generate  the  next  gain 
trajectory  plot  for  monitor  output.  It  is  important  to  note  that  the  gain  trajectories  are 
plotted  against  the  rea/  time  discrete  index,  k.  This  means  that  the  first  gains  calculated 
are  those  on  the  rightmost  edge  of  the  plot  while  the  first  gains  implemented  are  those 
on  the  leftmost  edge  of  the  plot.  Thus,  the  term  'steady  state'  as  it  applies  to  optimal 
feedbad^  gains  refers  to  the  zero-slope  property  of  the  left  side  of  the  gain  trajectory 
plot.  The  two  gain  trajectory  plots  for  the  example  problem  are  shown  in  Figures  A,3a 
and  A. 3b.  When  all  n  gain  plots  have  been  displayed  on  the  monitor  and/or  have  been 
printed  as  hardcopies,  the  program  continues  with  Screen  32. 

Screen  32 

Do  you  want  to  chongo  tho  NUfMR  OF  STAQCS  ? 

(AnoMor  y  or  n)  n 

If  the  user  is  not  satisfied  with  the  initial  choice  of  N,  then  a  new  value 
may  be  entered  at  this  time  by  answering  y  in  Screen  32.  OPTCON  presents  Screen  5 
for  this  purpose  and  subsequently  returns  to  the  sequence  beginning  with  Screen  28. 


DISCRETE  RERL  TIME  INDEX  Ik)  *10 


FEEDBACK  GRIN  (F2)  FOR  STATE  X2 

IMiniiizalion  over  RLL  STAGES), 


DISCRETE  REAL  TIME  INO 


The  most  likely  reason  for  the  user  to  take  advantage  of  the  option  in 
Screen  32  is  that  the  gain  trajectories  do  not  reach  steady  state  values  in  the  allotted 
number  of  time  intervals.  By  increasing  N,  the  user  may  be  able  to  force  the  gains  to 
reach  steady  state.  Since  the  gain  trajectories  in  Figures  A.3a  and  A.3b  demonstrate 
steady  state  properties,  there  is  no  need  to  change  N  in  Screen  32. 

6.  The  Time  Response 

Block  6  in  Figure  A.l  involves  calculation  of  the  time  response  of  the 
system  based  on  the  optimal  gains  computed  in  Block  3.  The  first  option  available  to 
the  user  in  this  section  is  the  phase  plane  graph  of  Xj  verses  X2.  Screens  33  through  36 
represent  the  program  sequence  that  results  when  a  phase  plane  is  requested  with  the 


following  constraints  : 

V 

10  seconds 

Xj(0) 

-  0  - 

Initial  condition  on  state  Xj 

XjfO) 

-  0  - 

Initial  condition  on  state  X2 

T(l) 

-  1  - 

Step  forcing  function  on  state  Xj 

r(2) 

-  0  - 

Ramp  forcing  function  on  state  x. 

Screen  33 


Do  you  Miin«  to  SM  •  PHASE  PLANE  of  XI  .vs.  X2  ? 
(AniMwr  y  or  n) 


Screen  34 

For  hoM  osny  sooonds  ?  10 
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Screen  35 

Ent«r 

tha  mlmm 

■nts  of  tho  INITIAL  STATE  vaetor  -  MUO) 

XKO) 

■  T 

0 

X2(0) 

■  ? 

0 

The  next  option  available  is  to  select  the  method  by  which  the  optimal 
gains  are  to  be  implemented.  Two  choices  are  available. 


Screen  37 

Salwt  •  0aln  knmmr  1  or  2  ) 

1)  Um  steady  state  0iin«  ovor  oil  o<«|>o  . 
Zi  Um  dynamic  eoino  . 

- ANSNER . 


If  the  first  option  is  chosen,  then  the  state  trajectories  are  calculated  using 
Equations  2.14  through  2.17  such  that  the  last  gain  matrix  calculated,  F(N),  is 
substituted  into  Equation  2.17  during  every  cycle  of  the  iteration  process.  The  user 
must  be  aware  of  this  procedure  when  selecting  option  1  in  Screen  37  because 
OPTCON  makes  no  attempt  to  verify  that  the  gains  have  indeed  reached  steady  state. 
If  the  user  selects  option  1  when  the  gains  do  not  exhibit  steady  state  properties,  then 
the  solution  is  not  optimal.  If  the  gain  trajectories  do  arrive  at  steady  state  prior  to  the 
N***  stage,  then  selection  of  option  1  in  Screen  37  may  be  appropriate.  The  time 
response  obtained  in  this  fashion  represents  the  behavior  of  the  system  using  a  fixed 
gain  feedback  scheme. 
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The  second  option  in  Screen  37  causes  the  feedback  gains  to  be 
dynamically  implemented  in  the  reverse  order  that  they  are  calculated.  The  user  is 
cautioned  that  such  implementation  may  not  yield  an  acceptable  time  response.  In  the 
example  problem,  the  gains  reach  steady  state  after  approximately  500  stages.  This 
corresponds  to  tf  -  5  seconds  when  At  -  0.01.  Consider  the  case  of  a  sampled  system 
which  has  a  transient  time  response  longer  than  S  seconds.  The  use  of  a  dynamic  gain 
schedule  would  be  disasterous  in  this  situation.  Because  the  gains  progress  towards 
zerp  as  t  approaches  5  seconds,  the  feedback  channel  is  gradually  ,  eliminated  from  the 
system.  The  slow  system,  however,  docs  not  have,  enough  time  to  reach  steady  state 
before  the  feedback  gains  go  to  zero.  The  error  signal  increases  without  bound  and  the 
system  rapidly  becomes  unstable.  Two  simple  solutions  for  such  a  situation  are  : 

1.  Increase  the  number  of  time  intervals,  N.  This  causes  the  steady  state  portion 
of  the  dynamic  gain  schedule  to  become  more  predominate. 

2.  Implemeni  steady  state  gains  instead  of  a  dynamic  gain  schedule. 

After  a  gain  schedule  is  selected  in  Screen  37,  OPTCON  begins  to  compute 
and  save  the  state  trajectories  for  x^  and  X2  using  Equations  2.14  through  2.17.  The 
message  in  Screen  38  informs  the  user  that  the  program  is  still  executing. 


Screen  38 

Caloulating  Plotting  Dot* 


After  the  state  trajectories  are  computed.  Screen  39  prompts  the  user. 


Screen  39 

READY  TO  OZSPUY  DRAHZNO 
Striko  any  Koy  to  eontinuo. 


The  monitor  is  cleared  upon  any  keystroke  and  the  Xj  verses  X2  phase  plane 
subsequently  appears  on  the  screen.  The  graph  remains  on  the  screen  until  the  user 
strikes  any  key.  The  monitor  then  clears  and  Screen  40  appears. 
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If  a  A7  is  entered  in  Screen  40,  then  the  program  advances  to  Screen  41. 
Otherwise,  the  message  in  Screen  38  reappears.  After  a  short  delay,  a  hardcopy  graph 
of  the  phase  plane  is  automatically  generated  on  the  printer.  See  Figure  A.4.  The 
program  then  advances  to  Screen  41. 

Screen  41 

Do  you  Mont  to  mo  a  tlM  nioporwo  of  your  ayotao  ? 

(AnsMor  y  or  n) 

If  a  is  entered  in  Screen  41,  then  the  first  run  of  OPTCON  is  complete. 
The  program  advances  to  Screen  44.  If  a>'  is  entered  in  Screen  41,  then  the  program 
prompts  the  user  to  enter  parameters  for  the  time  response.  Refer  to  Screens  34 
through  37.  It  is  not  required  that  the  user  enter  the  same  information  for  the  time 
response  that  was  entered  for  the  phase  plane.  OPTCON  recomputes  the  time 
response  on  every  run.  It  is  suggested,  however,  that  the  user  carefully  note  the 
parameters  that  are  entered  for  each  run.  Initial  conditions  and  command  inputs  are 
recorded  in  the  OPTFILE  but  this  information  does  not  appear  on  the  graphs.  After 
the  gain  schedule  is  selected  in  Screen  37,  OPTCON  begins  to  compute  and  save  the 
state  trajectories.  When  the  calculations  are  complete.  Screen  42  appears. 

Screen  42 

Do  yxi  mni  to  ms  toa  iia*  rmponM  tobla  on  to*  ooroan  7 
(Ansnar  y  or  n) 


X2  STATE 

-0.30  '  0.03  0.37  0.70  l.OU 


XI  vs.  X2  PHASE  PLRNE 


IMlniniiaHon  over  ALL  STAGES) 
USING  STERDT  STATE  GAIN  SCHEDULE 


Figure  A.4  Phase  Plane  for  Second  Order  Integrator  Example 
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A  n  response  causes  the  program  to  begin  generating  data  for  the  time 
response  plots.  The  user  is  cautioned  that  answering  y  in  Screen  42  may  result  in  a 
lengthy  delay  as  the  N  rows  of  data  are  scrolled  onto  the  monitor.  The  option  to  view 
this  data  on  the  monitor  exists  so  that  the  user  may  gather  exact  numerical  data 
without  exiting  OPTCON  to  examine  the  OPTFILE.  A  short  segment  of  the  tabular 
data  appears  in  Screen  43.  In  the  case  of  the  example  problem,  this  table  continues 
until  1000  rows  are  displayed. 


REAL 

TXM 

ZNoex 

REAL 

TXHE 

X(l) 

Screen  43 

Xf21 

X 

.0100 

.0000 

.0000 

2 

.0200 

.0000 

.0099 

s 

.0300 

.0002 

.0197 

4 

.0400 

.0004 

.0292 

S 

.OMO 

.0006 

.0364 

A 

.0400 

.0012 

.0479 

7 

.0700 

.0017 

.0670 

a 

.0600 

.0024 

.0469 

» 

.0900 

.0031 

.0744 

10 

.1000 

.0036 

.0632 

When  the  last  row  of  the  state  trajectory  table  appears  on  the  monitor,  or 
if  tabular  output  is  not  selected  in  Screen  42,  then  OPTCON  begins  to  generate  data 
for  the  state  trajectory  plots.  Each  state  is  plotted  verses  real  time.  During  the 
calculations,  the  messages  in  Screens  38  and  39  prompt  the  user.  State  Xj  is  plotted 
first  and  the  n^**  state  is  plotted  last.  The  user  may  examine  each  individual  graph  on 
the  monitor.  By  striking  any  key,  the  user  clears  the  graph  from'  the  screen  and  the 
message  in  Screen  40  reappears.  If  the  user  does  not  desire  a  hardcopy,  then  a  n 
response  allows  the  program  to  process  data  for  the  next  time  response  graph.  If  a 
hardcopy  is  desired,  then  a  y  is  entered  in  Screen  40  and  the  procedure  follows  exactly 
as  before.  See  typical  time  response  plots  in  Figures  A.5a  and  A.5b.  Afler  all  n  states 
are  plotted,  the  program  comnletes  Block  8  in  Figure  A.l.  The  first  run  of  OPTCON  is 
now  complete  and  the  user  must  answer  y  in  Screen  44  in  order  to  remain  in  the 
program.  If  a  n  is  entered  in  Screen  44,  then  execution  terminates  and  the  user  is 
immediately  returned  to  the  DOS  environment. 
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XI  TIHE  RESPONSE  BASED  ON 

IM^nialzalion  over  fiLL  STAGES) 
USING  OTNflHIC  CfllH  SCHEDULE 


AaoisHCdai  3iuis 


Figure  A.5a  State  Xj  Time  Response  for  Second  Order  Integrator 
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0.00  2.00  4.00  6.00  8.00 

REAL  TIHE  (sec) 


X2  TIHE  RESPONSE  BASED  ON 

(Hiniilz&i  Ion  over  RLL  STAGES) 


08'0  09'0  oro  02‘0  OO’O 

AU0133rbyi  3iblS 


Figure  A.5b  State  Xj  Time  Response  for  ^^cond  Order  Integrator 
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TIHE  (sec) 


*1^  ■ 

't’.v 


‘'pV? 


i>%  t", 


Assuming  that  the  user  desires  to  remain  in  OPTCON,  a  is  entered  in  Screen  44.  The 
next  section  of  the  program  is  referred  to  as  "the  main  menu"  and  is  demonstrated 
in  Screen  45. 


Screen  45 

SCLECT  ONE  OF  THE  FOLLQNZNB  OPTIONS! 

11  ChMiga  «ha  NUMBER  oi  STASES . 

..N 

2)  Chang*  th*  TERMINAL  atata  Mals^ting  aatrlx... 

..H 

3)  Changa  iha  TRAJECTORY  aiata  Maightlng  Matrix. 

..S 

41  Chang*  tha  CONTROL  tMightlng  factor. . 

S)  Changa  th*  praaant  A  and  S  aatrioa* 

..R 

41  Changa  tha  SAMPLE  INTERVAL . 

7)  Changa  th*  praaant  PHI  and  DEL  aatrioa* 

SI  Input  an  antiraly  NEN  SYSTEM 

91  NO  CHANSeS...RUN 

.DT 

101  EXIT  th*  prograa 

SELECTION. . .(  MUST  ba  a  nuabar  batmaan  1  and  10  1 

> 

It  is  not  necessary  to  describe  in  detail  the  operation  of  the  main  menu. 
The  user  should  enter  the  integer  value  that  applies  to  the  particular  modification  to  be 
made.  If  one  of  the  first  seven  options  is  selected,  the  program  responds  as  follows : 

1.  Echo  the  current  value(s)  of  the  parameter($)  to  be  modified. 

2.  Allow  the  user  to  keep  or  modify  the  selected  parameter(s). 

3.  Return  to  the  main  menu  for  further  modification,  continued  execution,  or 
termination  of  the  program. 

If  option  8  in  the  main  menu  is  selected,  then  OPTCON  returns  to  Screen  4 
and  allows  the  user  to  enter  new  information  for  all  parameters.  In  this  case,  no 
previous  values  are  remembered  by  the  program  and  execution  proceeds  just  as  if  this 
is  the  first  run.  The  OPTFILE,  however,  retains  all  information  from  any  previous 
runs. 


f  r* 
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If  option  9  is  selected  in  the  main  menu,  then  the  current  values  for  all 
system  and  cost  function  parameters  are  written  into  OPTFILE  to  signal  the  start  of  a 
new  run.  Program  execution  begins  with  the  gain  calculation  sequence  represented  by 
Blocks  in  Figure  A.  1.  Screen  28  applies.  The  user  may  rapidly  skip  through  the 
intermediate  steps  of  the  program  by  answering  n  to  several  consecutive  questions. 
For  instance,  suppose  that  the  user  changes  a  single  parameter  by  selecting  one  of  the 
first  seven  options  in  the  main  menu.  In  order  to  determine  the  effect  of  the  changed 
parameter  on  the  time  response  of  the  system,  the  following  sequence  of  messages  and 
responses  is  used. 

Screen  46 

SeieCT  QNi  OF  THE  FOLLOHXNe  OPTZQNSi 


1)  ChMiga  ^  NUHBER  of  STASES . N 

2)  Chang*  Ui*  TENHXNAL  ataia  tMtighting  Matrix . H 

3)  Changa  tha  TRAJECTORY  atata  malghting  Matrix... S 

A)  Changa  tha  COHTROL  waighting  factor . R 

5)  Changa  tha  praaant  A  and  B  aatrioaa 

6)  Changa  tha  SAHPLE  INTERVAL . OT 

7)  Changa  tha  praaant  PHI  and  DEL  aatrieas 

8)  Input  an  antiraly  NEN  SYSTEM 

9)  NO  CHANeES...RUN 
10 J  EXIT  tha  prograa 


SELECTION...!  MUST  ba  a  nuMbar  batuacn  1  and  10  ) . 9 

Do  you  want  to  saa  tha  gain  sohadula  tabla  on  tha  toraan  ? 
(Ail9i«ar  y  or  n)  n 

Oo  you  Mant  to  aaa  tha  gaina  plotted  7 

!AnaMar  y  or  n)  n 

Oo  you  want  to  aaa  a  PHASE  PLANE  of  XI  .va.  X2  7 
(Anawar  y  or  n)  n 

Oo  you  want  to  aaa  a  tiaa  raaponaa  of  your  aystao  7 
(Anawar  y  or  n)  y 


At  this  point,  the  user  may  examine  the  system  time  response  and  evaluate 
the  impact  of  the  newly  modified  parameter. 


By  selecting  option  10  in  the  main  menu,  the  user  is  allowed  to  exit  the 
program.  The  message  in  Screen  44  reappears  as  a  safety  mechanism  to  prevent 
inadvertant  ejection  from  the  program.  If  y  is  entered  in  Screen  44,  then  the  main 
menu  reappears  and  program  execution  continues.  Otherwise,  the  program  terminates 
and  control  is  returned  to  DOS. 

7.  OPTCON  Summary 

The  OPTCON  program  is  designed  specifically  so  that  the  user  can  easily 
modify  problem  parameters  and  rapidly  obtain  information  about  the  effects  of  those 
changes.*  Tabular  and  graphical  information  is  available  both  on  the  ihonitor  and  in 
hardcopy  form.  In  an  effort  to  make  the  program  user-friendly,  four  techniques  are 
employed  : 

1.  Menu  driven  options  prevail. 

2.  User  input  is  screened  for  valid  format. 

3.  User  inputs  are  echoed  on  the  monitor. 

4.  All  data  is  written  into  an  external  file  for  later  examination. 

The  OPTCON  program  is  quite  useful  as  an  interactive  design  tool  for  optimal  control 
systems.  Extensive  use  is  made  of  its  assets  during  the  design  of  an  optimal  controller 
for  the  AROD  in  chapter  three. 


133 


APPENDIX  B 

OPTCON  MAIN  PROGRAM  LISTING 


The  following  code  is  written  in  MICROSOFT  Fortran  and  is  intended  to  be 
used  on  an  IBM  compatible  system.  This  is  the  main  program  for  OPTCON  and  must 
be  linked  with  the  subroutines  found  in  Appendices  C  and  D.  In  addition,  the  Fortran, 
.Math,  and  PLOT88  libraries  must  be  linked  during  the  creation  of  an  executable  file. 


NOfloatcalls 
NO debug 


LL63.F0R 


LAST  MOD  12JULy87 


OK  SDL 

NEW  selective  state  plotting 
NEW  state  table  formatting 


COMMON  /BLKl/  A, B, PHI, DEL 
COMMON  /BLK2/  BEGTIM , FINTIM , NPTS , 

XNAML , YNAML , PNAMIL , PNAM2L , PNAM3L 
COMMON  /BLK3/  VTIME,VTIMSS,Vy,VYSS,VXXSS,VXYSS 
COMMON  / BLK4/  KFINAL , NSTAGE , NSTP 1 , ORDERN , GNSKED , USERGN , FNEG , 
INPUT , DT , AVG , AVG2 , MAXVAL , NINPTS 
KNAME , YNAME , PNAMEl , PNAME2 , PNAME3 
ION , ORDERN , IGOOD , CODE , NSTAGE , NSTPl , KFINAL , KPRIME , 


COMMON  .  _ 
INTEGER*2 


/BLK5/ 
*2  OPT] 


GNSKED , NPTS , lOPORT , MODEL , XNAML , YNAML , NCHARl , NCHAR2 , 
NCHAR3 , STVAR , I , J , SKIP , OK , SYSTEM , GAIN , DTFLAG , PLTYPE , 
CHNGN , SCREEN , NINPTS , NINPn , ORDNPl , GAINCH , GNSKD3 , 
STPLOT,PLOTCH 


,INPUT(8,1), 

FNEC(1006,8), 


VY( 1002) , VTIME  < 1002 ) , TFTEMP , TFINAL , DT , TIME , PNAMIL , PNAM2L , 
PNAM3L,WSS(9)  ,VTIMSS(9).VXXSS(9)  ,VXYSS(9)  ,AVG(8)  ,AVG2(8), 
MAXTIM , MAXVAL ( 6 ) , USERGN  6 , 8 ) 


CHARACTER*2  TEMP 
CHARACTER*3  ANS 
CHARACTER*20  NAME 
CHARACTER*30  XNAME. YNAME 
CHARACTER*51  PNAMEl ,PNAME2 ,PNAME 3 
CHARACTER*5  HDG(8‘ 

CHARACTER*4  HDG 
1 


HDG 
HDG 
HDG 
HDG 
HDG 
HDG 
HDG 
HDG  (8 
HDG, 
HDG2 
HDG2 
HDG2 
HDG2 
HDG2 
HDG2 
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^  ...^1  I  I  IS  ifLarU  tfVMU 


c 


HDG2(8)  »  'X(8)' 

0PEN<9 , FILE* ' OPTFILE ' , STATUS* ' NEW ) 


C 
C 

c****************  PRINT  OPTCON  HEADING  and  INPUT  PROBLEM  ID  ************ 
^*********************************************************************** 
C 

WRITE (*,2000) 

WRITE(*,2010) 

PAUSE 

WRITE (*,2015) 

READ  (*,2020,END*1S30)NAME 
C 

Q*********************************************************************** 

C*****^**********  HEADING  INFO  FOR  THE  OUTPUT  FILE  ■******■#***** 

C*A*************A**********)lt***3t*********************>lc:>*)lt********A****** 


c 

WRITE (9, 2030) 

WRITE (9, 2040) 

WRITE (9, 2050 )NAME 
WRITE (9, 2030) 

C 

C********************************************************************** 

c****************  ENTER  PLOTTER/PRINTER  MODEL  TYPE  ************ 
(]********************************************************************** 
C 

5  WRITE(*,2055) 

READ  (*,2070)TEMP 

CALL  COMPARE ( TEMP ,1,2, CODE , IGOOD ) 

IF(CODE.EQ.O)GOTO  5 
IF(IG00D  .EQ.  1)THEN 
lOPORT  »  0 
MODEL  >  1 
ELSE 

lOPORT  ■  0 
MODEL  «  60 
ENDIF 
C 

(j********************************************************************** 

C************  INITIALIZE  B  .DSL.USERGN  MATRICES  ************ 


DO 


6  I  » 
DO  6 


CONTINUE 


1.8 

J  “  1.8 
BU.J.'f 
DELh,J) 
USERGN(I, 


J) 


0.0 

0.0 

0.0 


C 

Q************************************************************'‘  ********* 

c****************  enter  THE  ORDER  OF  THE  SYSTEM  ************ 

C********************************************************************** 


C 

C****************  RESET  FINAL  , GNSKED , GAINCH , GNSKD3  ************ 

10  FINAL  ■  0 
GNSKED  -  1 
GAINCH  *  0 
GNSKD3  «  0 
WRITE (*,2060) 

READ  (*,2070) TEMP 

CALL  COMPARE ( TEMP ,1.8, CODE , IGOOD ) 

IF7cODE.EQ.O)OOTO  10 
ORDERN  «  IGOOD 
ORDNPl  ■  IGOOD  +  1 

C 

(^**************************************A******************************* 

c****************  enter  THE  NUMBER  OF  CONTROL  INPUTS  ************ 

C************************>f  ********************************************* 


c 


Ij  WRITE (*,2075) 

READ  (*,2070) TEMP 
CALL  COHPARE(TEMP,1.Q,CODE,IGOOD) 

IF(CODE.EQ.O)GOTO  15 
NIOTTS  ■  IGOOD 
NINPPl  a  IGOOD  +  1 
C 

C**^*************  ECHO  NINPTS  ************ 

C 

WRITE (*, 2076 )NINPTS 
C 

C****************  MODIFY  NINPTS  IF  NEEDED  ************ 

C 

16  WRITE(*,2077) 

READ  (*,2190) ANSWER 

IF(ANSWER.EQ.'N' .OR.ANSWER.EQ.'n')GOTO  17  ~- 

IF(ANSWER.EQ. 'Y' . OR. ANSWER. EQ. ' y' )GOTO  15 
GOTO  1 6 

17  CONTINUE 
C 

C******  SKIP  COST  FUNCTION  ENTRY  IF  NUMBER  OF  CONTROLS  .OT.  1  ******* 

C 

1F(NINPTS  .GT.  1)  THEN 
GNSKED  =  3 
GOTO  340 
ENDIF 

Q******************************************************A*************** 

C****************  ENTER  THE  NUMBER  OF  TIME  INTERVALS  ************ 

^***************************A***********************************A****** 

c 

20  WRITE (*,2080) 

READ  (*,*)NSTAGE 
IF(NSTAGE  .GT.  1000) GOTO  20 
NSTPl  a  NSTAGE  +  1 
IF(CHNGN  .EQ.  l)GOTO  780 
IF(FINAL  .EQ.  l)GOTO  1520 


^****************************A****************A******************A*A**** 
C**A***A**AA*AA*A  INPUT  THE  0  MATRIX  aaaaaaaaaaaa 
CA*AA**A*AAA*%*AA*AAAAAAA*AAAAA*AA*A**AAAAAA*A**A**AA*A**AA**AA*AA*AA**A 


30 -LOOP  »  0 

WRITE (*,2090) 

READ  (*,2070)TEMP 

CALL  COMPARE ( TEMP ,1,3, CODE , IGOOD ) 

IF ( CODE. EQ.O) GOTO  30 

OPTION  a  IGOOD 

GOTO (40, 50, 60)  OPTION 

GOTO  30 

40  WRITE (*,2100) 

GOTO  80 

50  WRITE(*,2110) 

GOTO  80 

60  WRITE(*,2120) 

80  DO  90  I  a  1,6rDERN 

DO  90  J  a  l,ORDERN 


IF  (I  .EQ^.  J)  THEN 
"T (OPTION  .EC 


IF 

IF  (OPTION  .E( 

IF (OPTION  .EL 

WRITZ(*, 2130)1, w 


I,J)  a  1.0 
I,J)  a  O.C 
“N 


ELSE 


READ 
ENDIF 


90  CONTINUE 


ENDI 


I,J)  a  0.0 


0A*AA********A**A 


ECHO  THE 


MATRIX 


AAAAAAAAAAAA 
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100  CONTINUE 

WR1TE(*,2140) 

DO  110  Z>l,ORDERN 


Q  MATRIX  IF  NEEDED 


MODIFY  THE 


MATRIX  AAAAAAAAAAAA 


qAAAAAAAAAAAAAAAA 
C 

120  WRITE (*,2160) 

READ  (*,2070)TEMP 
CALL  COMPARE (TEMP, 1.3, CODE, IGOOD) 

IF (CODE. EO.Oj GOTO  120 
OPTION  «  IGOOD 
GOTO (130 , so , 160 )OPTION 
GOTO  120 
C 

CAAAAAAAAAAAAAAA*  CHANGE  ONE  ELEMENT  OF  THE 
C 

130  WRITE(*,2170) 

READ  (><,*)I,J 

IF(I.LT.l  .OR.  I.GT.ORDERN  .OR.  J.LT.l  .OR.  J.GT.ORDERN)GOTO  130 

WR1TE(*, 2130)1, J 

READ 

write!* '2140) 

DO  140  Ial,ORDERN 

140  WRITE ( * , 2150 ) (Q ( I , J) , J*1 , ORDERN) 

150  WRITE (*,2180) 

READ  (^,2 190) ANSWER 

IF(ANSWER.EQ. 'N* .OR. ANSWER. EQ. 'n' )GOTO  160 
IF ( ANSWER. EQ. 'Y' .OR. ANSWER. EQ. 'y' )GOTO  130 
GOTO  ISO 

160  IF ( FINAL. EQ.l) GOTO  1520 


CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 
CAAAAAAAAAAAAAAAA  INPUT  THE  H  MATRIX  AAAAAAAAAAAA 

qAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

c 

170  LOOP  *  0 

WRITE (*,2200) 

READ  {*,2070)TEMP 

CALL  COMPARE ( TEMP ,1,3, CODE , IGOOD ) 

--IFXCODE.EJ.O)GOTO  170 


OPTION  »  IGOOD 
GOTO( 180, 190,200)  OPTION 
GOTO  170 

180  WRITE(*,2210) 

GOTO  210 

190  WRITE (*,2220) 

GOTO  210 

2C0  WRITE(*,2230) 

210  DO  220  I  «  1, ORDERN 

DO  220  J  -  1, ORDERN 

IF(I  ,e0.  J)  THEN 
IF  (OPTION  .E( 

IF (OPTION  .E( 
IF(OPTION  .E( 

WRITE (*,2*,w,*,« 
READ  (*,*)H(t,0) 
ENDIF 
ELSE 

0(1, J)  ■  0.0 
ENDIF 

220  CONTINUE 


I,J)  -  1. 
I,J)  »  0. 
N 


CAAAAAAAAAAAAAAAA 

C 

230  CONTINUE 

WRITE(*,2250) 

‘  ■  i»i,ofo 


ECHO  THE 


MATRIX 


AAAAAAAAAAAA 


DO  240 


ERN 
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240 


MODIFY  THE  H 


MATRIX  IF  NEEDED 


C 

250  WRITE (*,2160) 

READ  (><, 2070) TEMP 
CALL  COMPARE ( TEMP ,1,3, CODE , IGOOD ) 

IF7cODE.EQ.O)GOTO  250 
OPTION  ■  IGOOD 
GOTO (260 , 190 , 290 )0PT10N 
GOTO  250 
C 

C*******^********  CHANGE  ONE  ELEMENT  OF  THE 
C 

260  WRITE(*,2170) 

■‘•READ  (’*,*)!, J 

IF(I.LT.l  .OR.  I.GT.ORDERN  .OR.  J.LT.l  .OR.  J.GT.OROERN)GOTO  260 
WRITE (*, 2240 )1,J 

.i) 

TERN 


H  MATRIX  ************ 


WRITE(*,2250) 
DO  270  i-i.oto; 


270  WRITE(*,2150)(H(I,J),J«1,ORDERN) 


280  WRITE (*,2180) 

READ  (^,2 190) ANSWER 

IF(ANSWER.EO. 'N* .OR. ANSWER. E 
IF ( ANSWER. EQ. 'Y* .OR. ANSWER 
GOTO  280 

290  IF ( FINAL. EQ.1)GOTO  1520 


.EQ.'n')OOTO  290 
.EQ.'y'jGOTO  260 


C*********************************************************************** 

^***A************  INPUT  R  ************ 
^*********************************************************************** 
C 


300  WRITE (*,2260) 
READ  (’<,*)R 


ECHO 


^**************** 

C 

310  WRITE (*, 2270 )R 
C 

(j**************** 

C  -  — 

320  WRITE (*,2280) 

READ  (’^, 2 190) ANSWER 

IF(ANSWER.Eg. 'N* .OR. ANSWER. EQ. 'n'}GOTO  330 
l.EQ.'Y'.OR.ANSWER.EQ.'' 


MODIFY 


IF  NEEDED 


************ 


************ 


IF (ANSWER 
GOTO  320 
330  IF ( FINAL. EQ.l) GOTO  1520 


'y')GOTO  300 


C*********************************************************************** 

Q****************  CHOOSE  TO  ENTER  EITHER  A  ************ 

C****************  CONTINUOUS  TIME  SYSTEM  OR  A  ************ 

C****************  DISCRETE  TIME  SYSTEM  ************ 

Q*********************************************************************** 

C 

340  WRITE (*,2290) 

READ  (^,2070)TEMP 

CALL  COMPARE (TEMP, 0,1, CODE, IGOOD) 

IF(CODE.Eg.O)GOTO  340 
SYSTEM  ■  IC 


[GOOD 


350 


360 


IF(SYSTEM)350 , 350 , 360 
WRITE (*,2300) 

READ  (*,2190)ANSWER 

IF(ANSWER.EQ. 'Y* . OR. ANSWER. EQ. ' y* )GOTO  370 
SYSTEM  ■  1 
WRITE (*,23 10) 

READ  (*,2190)ANSWER 

IF ( ANSWER. EQ. 'Y' .OR.ANSWER.EQ. 'y* )GOTO  590 
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SYSTEM  >■ 
GOTO  350 


(;A>.'AAAAA*AAAAi<(A*ili  INPUT  THE  A  MATRIX  ************ 

Q********************************************************************** 

C 

370  WRITE(*,2320) 

DO  380  I«l,OROERN 
DO  380  J»l,ORDERN 
WRITE (*,2330) I, J 
READ  (*,*)A(I,.J) 

380  CONTINUE 
C 

C****************  DO  NOT  ALLOW  CHANOKS  TO  A  «nd  B  ************ 
C******j^*********  IF  A  DISCRETE  TIME  SYSTEM  WAS  ENTERED  ************ 

390  CONTINUE 

IF (SYSTEM  .EQ.  1)THEN 
WRITE(*,2333) 

READ  (*,2070)TEMP 
GOTO  1520 
ENDIF 


0****************  ECHO  THE 

C 

WRITE (*,2340) 

DO  400  I«l,0toERN 

400  WRITE(*,2150)<A(I,J),J*1,ORDERN) 
C 

0****************  MODIFY  THE  A 

C 


MATRIX 


MATRIX  IF  NEEDED 


************ 


************ 


410  WR1TE(*,2160) 

READ  (^,207(5) TEMP 


CALL  COMPARE ( TEMP ,1.3, CODE , IGOOD ) 
IF7cODE.EQ.O)GOTO  410 
OPTION  -  IGOOD 
GOTO (420 , 370 , 450 )0PT10N 
GOTO  410 


CHANGE  ONE  ELEMENT  OF  THE 


MATRIX  ************ 


C**************** 

C 

420  1«ITE(*,2170) 

READ 

IF(I.LT.l  .OR.  I.GT.ORDERN  .OR.  J.LT.l  .OR.  J.GT.ORDERN)GOTO  420 
WRITE?*, 2330) I, J 

--  p, 

}ERN 


READ 
WRITE^ 
DO  430 


: 


:2i4dr 

i»i,otei 

-njr- 


430  WRITE (*,2150)(A(I,J),J«1, ORDERN) 
440  WRITE (*,2180) 

READ  (*,21 90) ANSWER 
IF(ANSWER.EQ. 'N* . 


SO. 

SQ. 


IF ( ANSWER. E 
GOTO  440 
450  IF(FINAL.EQ.l)GOTO  480 


.OR. ANSWER. E 


Y' .OR. ANSWER 


.EO.'n' 
.EQ. 'y' 


)GOTO  450 
)GOTO  420 


C*********************************************************************** 

C****************  INPUT  THE  B  MATRIX  ************ 

^********************************A***************A*******************A** 

C 

460  WRITE (*,2350) 

DO  470  I«l, ORDERN 
DO  470  J  ■  1,NINPTS 
WRITE (*, 2360)1, J 
READ  (*  *)B(I,jr) 

470  CONTINUE 
C 

Q****************  ECHO  THE  B  MATRIX  ************ 

C 
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480  CONTimJS 

WRITER*. 2370) 

DO  490  I»1,QRDERN 

490  WRITE(*, 2150) (8(1. J),J-1,N1MPTS) 

MODIFY  THE  B 


MATRIX  IF  NEEDED 


it****itAit***-k 


C****ilt****<»)(^****A 

C 

500  WRITE (*,2160) 

READ  (\2070)TEMP 
CALL  COMP ARE (TEMP. 1,3. CODE, IGOOD) 

IFTcODE.EQ.OjOOTO  500 
OPTION  «  lOOOD 
GOTO (5 1 0 , 460 , 540 ) OPTION 
GOTO  500 
C 

C*AA***A*>kiitAA*A*A  CHANGE  ONE  ELEMENT  OF  THE 
C 

510  WRITE(*  2170) 

READ 

IF(I.LT.l  .OR.  I.GT.ORDERN  .OR.  J.LT.l  .OR.  J.GT.NINPTS)GOTO  510 
IITe7*, 2360)1, J 


B  MATRIX  AAAAAAil^AAAAA 


WRITE (*,2360) I, J 
READ  (*  *)S(f,i) 
WRITE(*,2370) 

DO  520  I-l .ORDER 


DERN 

520  WRITE(*,2150)(B(I.J),J«1,NINPTS) 

530  WRITE (*,2180) 

READ  (^,2 190) ANSWER 

IF(ANSWER.EQ. 'N* .OR. ANSWER. EQ. 'n* )GOTO  540 
IF(ANSWER.EQ. 'Y' .OR. ANSWER. EQ.*'  '  * 

GOTO  530 

540  IF ( FINAL. EQ.DGOTO  1520 


'y'jGOTO  510 


^AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 
CAaaaaaaaaaaaaaaa  INPUT  THE  SAMPLE  TIMS....DT  aaaaaaaaaaaa 

^AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA^AAAymAAAMKAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

C 

550  WRITE (*,2380) 

READ  (*,*)DT 
DTFLAG  ■  1 
C 

CAaaaaaaaaaaaaaaa  If  A  DISCRETE  TIME  SYSTEM  WAS  ENTERED 
CAAAAAAAAAAAAAAAA  AND  NO  VALUE  FOR  DT  HAS  BEEN  ENTERED 
CAAAAAAAAAAAAAAAA  THEN  PRINT  OUT  A  MESSAGE 

C 

560  IF (DTFLAG  ,EQ.  0)THEN 
WRITE (*.2335) 

READ  (*,2070)TEMP 
GOTO  1520 
ENDIF  . 


AAAAAAAAAAA 

AAAAAAAAAAA 

AAAAAAAAAAA 


ECHO  THE  SAMPLE  TIME....DT 


aaaaaaaaaaaa 


AAAAAAAAAAAA 


CAAAAAAAAAAAAAAAA 

C 

WRITE (*, 2390 )DT 

CAAAAAAAAAAAAAAAA  MODIFY  THE  SAMPLE  TIME  IF  NEEDED 
C 

570  WRITE (*.2400) 

READ  (^,2190) ANSWER 

IF(ANSWER.EQ.'N'. OR. ANSWER. EQ.'n')GOTO  580 
IF(ANSWER.EQ. 'Y' . OR. ANSWER. EQ. 'y* )OOTO  550 
GOTO  570 
C 

CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 
CAAAAAAAAAAAAAAAA  CONVERT  A  and  B  TO  PHI  and  DEL  AAAAAAAAAAA 

CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

c 

580  IF (SYSTEM  .EQ.  0)  THEN 

CALL  PHIDEL(DT,0R0EP31.NINPTS) 

ENDIF 

IF(FINAL.EQ.1)GOTO  1520 
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jiBWiR  anil  mil 


GOTO  760 


C******A*********  INPUT  THE  PHI  MATRIX  *********** 

^*********************************************************************** 
C 

590  CONTINUE 
DTFLAG  *0 
600  WRITER, 2410) 

DO  610  I-l,ORDERN 
DO  610  J«l,OROERN 
WRITE(*, 2420)1, J 
READ  (*,*)PHl(i,J) 

610  CONTINUE 
C 

C****************  DO  NOT  ALLOW  CHANGES  TO  PHI  and  DEL  ************ 
C*****if**********  IF  A  CONTINUOUS  TIME  SYSTEM  WAS  ENTERED  ************ 
C 

620  CONTINUE 

^  0)THEN 
r42S) 


IF{SYSTEM^^.E^^. 


ENDIF 


WRITE{*, _ . 

READ  (*.2070)TEMP 
GOTO  1520 


ECHO  THE  PHI  MATRIX 


************ 


************ 


************ 


^A*************** 

c 

WRITE (*2430) 

DO  630  fal.ORDERN 

630  WRITE(*,2150)(PHI(I,J),J«1,ORDERM) 

C****************  MODIFY  THE  PHI  MATRIX  IF  NEEDED 
C 

640  WRITE (*2160) 

READ  (*,2070)TEMP 
CALL  COMPARE(TEMP, 1.3, CODE, IGOOD) 

IFTcODE.EQ.oIgOTO  640 
OPTION  ■  IGOOD 
GOTO (650,600,680) OPTION 
GOTO  640 
C 

C****************  CHANGE  ONE  ELEMENT  OF  THE  PHI  MATRIX 
C  — 

650  WRITE (*.2170) 

READ  (^,*)I,J 

IF(I.LT.l  .OR.  I.GT.ORDERN  .OR.  J.LT.l  .OR.  J.GT.ORDERN)GOTO  650 
WRITE?*, 2420)1, J 
READ  (*,*)PHI(I,J) 

WRITE (*,2430) 

00,660  I«l,OtoERN 

660  WRITE(*,2150}(PHI(I,J),J>1,OROERN) 

670  WRITE  *; 2180) 

READ  (*,2190) ANSWER 

IF(XnSWER.EQ. 'N* .OR. ANSWER. EQ. 'n* )GOTO  680 
IF( ANSWER. EQ. 'Y' . OR. ANSWER. EQ. )GOTO  650 
GOTO  670 

680  IF(FINAL.EQ.l)GOTO  710 
C 

Q*********************************************************************** 

C****************  INPUT  THE  DEL  MATRIX  ************ 

C*********************************************************************** 

C 

690  WRITE (*,2440) 

DO  700  I»1,OROERN 
DO  700  J«1,NINPTS 
WRITE(*,2450)I,J 
READ  (*,*)DELU,J> 

700  CONTINUE 
C 

C****************  ECHO  THE  DEL  MATRIX  ************ 
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710  CONTINUE 

WR1TE(*.2460) 

DO  720  I«l,ORDERN 

720  WR1TE(*,2150)(DEL(I,J).  J*1,N1NPTS) 

C 

C****************  MODIFY  THE  DEL  MATRIX  IF  NEEDED  *********)iti»t* 
C 

730  WRITE (*.2160) 

READ  (*,2070) TEMP 

CALL  COMPARE(TEMP,1.3,CODE,IOOOD) 

IFTcODE.EQ.OJGOTO  730 
OPTION  ■  IGOOD 
GOTO (740 , 690 , 770 )OPTION 
GOTO  730 

C  „ 

C*****A**********  CHANGE  ONE  ELEMENT  OF  THE  DEL  MATRIX  ************ 
C 

740  WRITE (*2170) 

READ  (*,*)!, J 

IF(I.LT.l  .OR.  I.GT.ORDERN  .OR.  J.LT.l  .OR.  J.GT.NINPTS)GOTO  740 
WRITE?*, 2450 )1,J 
READ  (*,*)DEL(1,J) 

WRITE (*,2460) 

DO  750  I«l,ORDERN 


750  WRITE(*,2150)(DEL(I,J),J«1,NINPTS) 
760  WRITE  *; 2180) 

READ  (^,2190) ANSWER 

IF(ANSWER.EQ. 'N* .OR.ANSWER.EQ. ' n* 
IF(ANSWER.EQ. *Y' .OR.ANSWER.EQ. 'y* 
GOTO  760 

770  IF(FINAL.EQ.l)GOTO  1520 


)GOTO  770 
JGOTO  740 


CAAAAAAAAAAAAAAAA  WRITE  ALL  CURRENT  INFORMATION  TO  THE  AAAAAAAAAAAA 
CAAAAAAAAAAAAAAAA  OUTPUT  FILE  AAAAAAAAAAAA 

^AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

c 

■  780  FINAL  ■  0 

WRITE (9, 2030] 

WRITE(9,2470)ORDERN 

- [9,2475)NINPTS 

3)  GOTO  805 
CAC 


WRITE 
-IF(GNSI«D  .EC 


,J),J>I,ORDERN) 


WRITE ( 9 , 2480  TNSTAGE 
WRITE (9, 2140] 

TRACEQ  -  0.0 
DO  790  Ial,ORDERN 

TRACEQ  «  TRACEQ  +  Q(I,I) 

790  WRITE (f, 2490 )(Q a, '  '  * 

WRITE(9,2i50) 

DO  800  I>l,0toERN 

800  WRITEj9,2490)(H(I,J),J»l,0RDERN) 
WRITE (9,2270)R 
805  IF(SySTEM)  810,810,840 
810  WRITE (9, 2340) 

DO  820  I«l,OROERN 

820  WRITE  (9 ,2490 )  ( A  ( I ,  J  ) ,  J»l.,  ORDERN) 

WRITE (§^2370^ 


830 


DO  830  !• 


*  ,  WnVM 


;DERN 


I,J),J»1,NINPTS) 


WRITE(9,249C , . 
write] 9, 2390 
840  WRITE (9, 2430) 

DO  850  I>1, ORDERN 

850  WRITE (9 , 2490 ) (PHI ( I , J) , J«1 , ORDERN) 

WRITE]9,2460) 

DO  860  fsl, ORDERN 

860  WRITE  (9.2490)(DEL(I,J),J*l,NINPTS) 
WRITE(9,2d30) 

IF(GNSKED  .EO.  3)  THEN 

CAAAAAaaa  no  OPTIMAL  GAINS  ARE  TO  BE  CALCULATED 


AAAAAAAAAAAA 
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GOTO  1010 
ENDIF 

IFTtRACEQ) 870 .870 , 880 
870  WRITE(9,«00) 

PNAHE2"  •(Klnimum  TERMINAL  STATES  Control)' 
PNAM2L*  33 
GO  TO  890 
680  WRZTE(9,2S10) 


PNAME2-  ' (Mlnlmiz«tlon  ovor  ALL  STAGES 
PNAM2L-  30 


) 


AAAAAAAAAAAA 


CAAAAAAAAAAAAAAA*  INITIALIZE  MATRICES  PRIOR  TO  AAAAAAAAA*** 

C*a*a***a*a*a**a*  calculating  OPTIMAL  GAINS  aaaaaaaaaaa* 

^AAAAAAAAAAAAAAAAAAAAAAAAAAAi^AAkAAAAAAAAAAAAkAAAAAAAAAAAAAAAAAAAAAAAAAA 

c 

890  Continue 

00  900  I«1.0R0ERN 
EM(I}  «  0.0 
FMU)  ■  0.0 
DO  900  J«l,OROERN 
GM(I,J)  »  0.0 
HM(I,J)  *  0.0 
900  P(I,J)-H(I,J) 

CAAAAAAAAA  QQ  YOU  WANT  TO  SEE  THE  GAINS  TABLE  ON  THE  SCREEN  ?  ********* 
C 

WRITE(*,2S15) 

READ  (*,2190) ANSWER 

.  IF(ANSWER.EQ.'N' .OR. ANSWER. EQ. •n')SCREEN  ■  0 
IF(ANSWER.EQ.'Y'.OR.ANSWER.EQ.'y')SCREEN  -  1 

C****************  print  HEADING  FOR  OUTPUT  TABLE 
C****************  OPTIMAL  GAINS 

C 

IF(SCREEN  .EQ.  1)THEN 

WRITER*, 2520) (HDG(I) ,I-l,ORDERN) 

WRITE (*,2030) 

ENDIF 

WRITE(9,2S20)(HDG(I),I«1,OROERN) 

WRITE (9, 2030) 

^********************************************************************** 

C****************  LOOP  TO  ITERATE  THE  RICATI  EQUATIONS  ********** 
^A*************************************^****************************** 

c 

DO  1000  KKal,NSTAGE 
KREAL  >  NSTPl  •  KK 
DEN>0.0 

DO  910  Ial,OROERN 
DO  910  J«1,0RDERN 
EM(I)  ■  EM(I)  +  DEL(J,1)  *  P(J,I) 

DO  930  I«x,OROERN 
DO  920  J«1,0RDERN 
FM(I)  -  FM(I)  +  EM<'J)  *  PHI(J,I) 

DEN  ■  DEN  +  EM(I)  *  DELU,!) 

DEN  >  DEN  -r  R 


910 


920 

930 


C*********  ensure  THAT  THE  DENOMINATOR  DOES  NOT  GO  TO  ZERO  *********** 

C 

IF(  DEN  .EQ.  0  )THEN 
WRITE  C*. 2530 )KK-1 
WRITE (9, 2530 )KK-1 
NSTAGE  ■  KK  -  1 
GOTO  1007 
ENDIF 
C 

C****************  CALCULATE  OPTIMAL  GAINS  FOR  THIS  STEP  ********** 
C 

DO  940  Ial,ORDERN 
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EM(I) 


FTRAN(I)  ■  -FM(X)/DEN 
FNBO<ia(,Z)  ■  FT^(Z) 
FM(1)  ■  0.0 


CAAAAAAAAAAAAiInliiltA  PRINT  OPTIMAL  OAIHS  FOR  THIS  STEP  ***i^***A**:fc 

c 

IF(SCREEN  .EQ.  IJTHEN 

IFCORDE^  .GT.  4)  THEN 

WRITE  ( ,  2540 )  KX ,  KREAL ,  ( FTRAN ( I ) ,  1*1 ,  ORDERN) 

EXiSE 

WRITE(*, 2S41)KK, KREAL, (FTRAN(I) ,1*1, ORDERN) 

SNDZF 

ENDIF 

IF<ORDERN  .OT.  4)  THEN 

WRITE ( 9 , 2 540 ) KK , KREAL , ( FTRAN ( I ) , 1*1 , ORDERN ) 

*  -  ELSE  “ 

WRITE ( 9 , 2541 ) KX , KREAL ,( FTRAN (I ), 1*1 , ORDERN) 

ENDIF 

C 

CALCULATE  PSI(K,I,J)  ****:******* 

^  DO  950  1*1, ORDERN 

DO  950  J*1 .ORDERN 

950  PSI(I,J)  *  PHI(I,J)  +  DEL(I,1)  *  FTRAN(J) 

C 

Q****************  CALCULATE  P  (K,I,J)  **:fc*i*f**A**** 

^  DO  960  1*1, ORDERN 

DO  960  Jal, ORDERN 
DO  960  L*l, ORDERN 

DO  980  J*l, ORDERN 
DO  970  L«l, ORDERN 

970  HM(I,J)  *  HMU/J)  +  GM(I,L)  *  PS1(L,J) 

«««  •  mi.J)  +  Qa,J)  +  R  *  FTRAN(I)  *  FTRAN(J) 

980  HM(1,J)  *  0,0 

DO  990  1*1, ORDERN 
DO  990  J*l, ORDERN 
990  GM(1,J)  *  0.0 

C****************  DISCRETE  TIME  VECTOR  FOR  PLOTTING  GAINS  ************ 


c 

VTIME(KK)  ■  KK 

1000  CONTINUE 
C 

C***************  DO  YOU  WANT  TO  SEE  THE  GAINS  PLOTTED  ?  *********** 

C 

1001  WRITE(*,254S) 

READ  (*,2190 'ANSWER 

IF ( ANSWER. EQ. 'N' .OR.ANSWER.EQ. 'n' )GOTO  1006 
IF(ANSWER.EQ. 'Y* .OR.ANSWER.EQ. 'y' )GOTO  1002 
GOTO  1001 
C 

C****************  LOOP  TO  PLOT  OUT  THE  GAINS  ************ 

C 

1002  DO  1005  GAIN  *  1, ORDERN 
C 

C****************  SET  THE  GAIN  PLOT  TITLE  ************ 

C 

IF  (GAIN.  E< 

IF  (GAIN.  E< 

IF(GAIN.E< 

IF(GAIN.Ei 
IF(GAIN.Ei 
IF  (GAIN.  El 
IF  (GAIN.  El 
IF(GAIN.Ei 


.1)PNAME1  *  'FEEDBACK  GAIN  (FI)  FOR  STATE  XI 
.2)PNAHE1  *  'FEEDBACK  GAIN  (F2)  FOR  STATE  X2 
.3)PNAHE1  *  'FEEDBACK  GAIN  (F3)  FOR  STATE  X3 
.4)PNAME1  *  'FEEDBACK  GAIN  (F4  i  FOR  STATE  X4 
.5)PNAME1  *  'FEEDBACK  GAIN  iF5)  FOR  STATE  XS 
.6)PNAME1  *  'FEEDBACK  GAIN  (F6)  FOR  STATE  X6 
.7)PNAME1  a  'FEEDBACK  GAIN  (F7)  FOR  STATE  X7 
.8)PNAHE1  a  'FEEDBACK  GAIN  (F8)  FOR  STATE  X8 
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PNAMIL  «  31. 

C 

CAAAAAAAAAAAAAAAA  PLOT  88  GRAPHICS  AAAAAAAAAAA 

qAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAkAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

c 

CAAAAAAAAAAAAAAAA  SET  UP  INITIAL  PARAMETERS  FOR  GAIN  PLOT  *aaaaaaaaaa 
BEGTIM  -  0.0 
FINTIM  »  NSTAGE 
NPTS  »  NSTAGE 

DO  1003  J  ■  1,7 
VYSS(J)  ■  0.0 

1003  VT1MSS(J)  ■  ((FINTIM  -  BEGTIM)/6. )*(J-1) 

CAAAAAAAAAAAAAAAA  GENERATE  GAIN  VECTOR  FOR  PLOTTING  GAINS  aaaaaaaaaaaa 
C 

*•  DO  1004  KREAL  ■  1, NSTAGE  ' 

KK  -  NSTPl  -  KREAL 

VY (KREAL)  «  FNEG?KK,GAIN) 

C*********  TEST  LINE  FOR  SELECTING  PROPER  COLUMN  OF  GAINS  FOLLOWS  **** 
CAAAAAAAAA  SEE  LL51  FOR  COMPILED  VERSION  **** 

C  WRITE (*,*)  GAIN, KREAL, KK,FNEO(KK, GAIN) 

1004  CONTINUE 
C 

CAAAAAAAAAAAAAAAA  MAKE  THE  GAIN  PLOT  AAAAAAAAAAA 

c 

IF (GAIN  .EQ.  1) PAUSE 
CALL  GRAPH(99,S9,1) 

CAAAAAAAAAAAAAAAA  JS  A  HARDCOPY  OF  THE  GAIN  PLOT  DESIRED  ?  AAAAAAAAAAAA 

c 

WRITE (*,2595) 

READ  (*,2190) ANSWER 

IF(ANSWER.EQ. 'Y' .OR.ANSWER.EQ. 'y' ) 

+  CJ^L  GRAPH(I0P0RT,M0DEL,1) 

CONTINUE 

1005  CONTINUE 
C 

CAAAAAAAAAAAAAAAA  DO  YOU  WANT  TO  CHANGE  NSTAGE  ?  aaaaaaaaaaaa 

c 

1006  CHNGN  >  0 

WRITE (*,2546) 

--READ  (*,21 90) ANSWER 

IF(ANSWER.EQ. 'Y' .OR. ANSWER. EQ. 'y* )THEN 
CHNGN  ■  1 
GOTO  20 

ELSE IF ( ANSWER. EQ. 'N' .OR.ANSWER.EQ. 'n* )THEN 
GOTO  1007 
ELSE 

GOT01006 

ENDIF 

C 

CAAAAAAAAAAAAAAAA  IS  ^  PHASE  PLANE  DESIRED  ?  AAAAAAAAAAAA 

C 

1007  WRITE(*,2547) 

READ  (*,2190) ANSWER 

IF(ANSWER.£Q.'Y'. OR. ANSWER. SQ.'y')  THEN 
NSTPl  ■  NSTAGE  1 
PHASE  »  1 
GOTO  1025 
ENDIF 

IF(ANSWER.EQ. 'N' .OR.ANSWER.EQ. 'n* )GOTO  1010 
GOTO  1007 
C 

CAAAAAAAAAAAAAAAA  IS  A  TIME  RESPONSE  DESIRED  ?  AAAAAAAAAAAA 

C 

1010  PHASE  «  0 

WRITE(*,2550) 

READ  (*,2190)ANSWER 

IF()(NSWER.EQ.  'Y'  .OR.ANSWER.EQ. '  y*  )GOTO  1020 
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I? (ANSWER. EQ. 'N* .OR. ANSWER. EQ. 'n* )OOTO  ISIO 
GOTO  1010 

Ci>tAAAAiitAiitiii«icA*AA**  qrj^H  IS  TO  BE  A  TIME  RESPONSE 

1020  NSTPl  -  NSTAGE  1 
PLTYPE  ■  3 
C 

C*AAA****A**AAA**  HQW  MANY  SECONDS  ? 

C 

1025  WRITE(*  2560) 

BEAD  (*,*)Tn 

C*AAA********A***  INPOT  DT  IF  NOT  ALREADY  KNOWN 

C 

IF(DTFLAG  .EO.  0)  THEN 
WRITE (*,2380) 

READ  (*,*)DT 
DTFLAG  ■  1 
ENDIF 
C 

CAAAAAAAAAAAAAAAA  CALCULATE  FINAL  VALUE  OF  K 

C 

KFINAL  »  NINT(TF1NAL/DT) 

TFTEMP  «  KFINAL  *  DT 
IF(TFTEMP  .LT.  TFINAL)  THEN 
KFINAL  ■  KFINAL  +  1 
TFINAL  ■  KFINAL  *  DT 
ENDIF 
C 

CAAAAAAAAAAAAAAAA  ENSURE  THAT  ENOUGH  GAINS  ARE  CALCULATED 
CAA*A*******A****  XO  COVER  THE  DESIRED  TIME  RANGE 

C 

IF  (GNSKED  .EQ.  3)  GOTO  1029 
IF((KFINAL-1)  .GT.  NSTAGE)  THEN 
MAXTIM  ■  DT  *  NSTAGE 
WRITE(*.2S61)MAXTIM 
GOTO  1025 
ENDIF 
C 

CAAAAAAAAAAAAAAAA  READ  IN  THE  INITIAL  STATE  VECTOR 
C 

1029  WR1TE(*,2565) 

DO  1030  lal.ORDERN 
WRITE(*, 2566)1 
READ  (*,*)Xk6(I,1) 

1030  CONTINUE 
C 

CA*A*****AAA*AAAA  READ  IN  THE  COMMAND  INPUT  VECTOR 
C 

WRITE(*,2570) 

DO  1035  I*l,OROERN 
WRITE(^,2580JI 
READ  (*,*)INPUT(I,1) 

1035  CONTINUE 
C 

C*********  WRITE  INITIAL  STATE  AND  COMMAND  INPUT  VECTOR 
C*********  XO  OUTPUT  FILE 

C 

WRITE (9, 2030) 

WRITE (9, 2584) 

DO  1036  I  »  i.ORDERN 

I,XK0(I,1),INPUT(I,1) 

1036  CONTINUE 


C*A********:*i*(***iHt 

Q*)lt***j*ti*:A:<t*AA*AA* 

C 

1040  IF(NINPTS  .EQ. 

IF(GAINa{ 


CHOOSE  EITHER  STEADY  STATE  GAINS  (1) 
OR  DYNAMIC  GAINS  (2) 

OR  USER  DEFINED  GAINS  (3) 

IF  ONLY  ONE  CONTROL  INPUT  IS  USED 


1)  THEN 
.NE. 


2)THEN 
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itiUtiiit****** 


************ 


********** 

********** 


************ 

************ 

************ 

************ 


AAAAAAi(%AAA'A 


AAAAAAAAAAAA 


WRITE (*,2590) 

READ  (*,2070)TEMP 
CALL  COMPARE (TEMP, 1.3. CODE, IGOOD) 

IFTcODE.EQ.OJGOTO  1040 
GNSKED  «  IGOOD 
ELSE 

GAINCH  «  1 
ENDIF 
ELSE 

GNSKED  a  3 
ENDIF 

GOTO (141. 142, 143)  GNSKED 

CAAAAAAAaaaaaaaaa  use  steady  STATS  GAINS 

141  PNAME3  a  'OPTIMUM  STEADY  STATE  GAIN  SCHEDULE' 

PNAM3L  a  34. 
goto  1054 

CAAAAAAAAAAAAAAAA  USE  DYNAMIC  GAINS 

142  PNAME3  «  'OPTIMUM  DYNAMIC  GAIN  SCHEDULE' 

PNAM3L  a  29. 

GOTO  1054 

CAAAAAAAAAAAAAAAA  IMPLEMENT  USER  DEFINED  FEEDBACK  GAINS  AAAaaaaaaaaa 

143  PNAME2  a  'Implementing' 

PNAM2L  a  12. 

PNAME3  a  'USER  DEFINED  GAINS' 

PNAH3L  a  13. 

IF(  GNSKD3  .EQ.  1  )  GOTO  1043 

i  <3nskd3  .eq.  i  )  goto  1043 

CAAAAAAAAAAAAAAAA  INPUT  USER  DEFINED  FEEDBACK  GAINS 

1044  DO  1045  I  a  1 .NINPTS 

DO  1045  J  a  l.ORDERN 
WRITE(*,2592)  I,J 
READ  (*,*)  ukRGN(I,J) 

1045  CONTINUE 

GNSK03  a  1 

CAAAAAAAAAAAAAAAA  ECHO  USER  DEFINED  FEEDBACK  MATRIX 

1043  WRITE(*.2593) 

DO  1046  lal. NINPTS 

1046  WRITE(*,2594)(USERGN(1,J),  Jal,ORDBRN) 

C 

CAAAAAAAAAAAAAAAA  MODIFY  THE  USER  DEFINED  GAINS  IF  NEEDED 

1047 "HRITE(*.2160) 

READ  (*,2070) TEMP 

CALL  COMPARE (TEMP, 1.3, CODE, IGOOD) 

I fTcODE . EQ . 0 ) GOTO  1047 
OPTION  a  looOD 
GOTO (1048 , 1044 , 1052 ) OPTION 
GOTO  1047 
C 

CAAAAAAAAA  CHANGE  ONE  ELEMENT  OF  USER  DEFINED  GAIN  MATRIX  **aaaaaaaaaa 

1048  WRITE(*,2170) 

REM  (*,*)I.J 

IF(I.LT.l  .OR.  I, GT. NINPTS  .OR.  J.LT.l 
+  _  .OR.  J.GT.ORDERN)GOTO  1048 

WRITE (*,2592) I, J 
READ  (*,*)USERGN(I,. 

WRITE (*,2593T 
DO  1049  lal, NINPTS 

1049  WRITE ( *  2 594 ) ( USERGN ( I , J ) , Jal , ORDERN ) 

1051  WRITE(*,2180) 

READ  (A, 2 190) ANSWER 


AAAAAAAAAAAA 


AAAAAAAAA 


-J) 


IF(ANSWER.EQ. *N' .OR. ANSWER. EQ. 'n 
IF(ANSWER.EQ. 'Y* .OR. ANSWER. EQ. 'y 
GOTO  1051 


')GOTO  1052 
')OOTO  1048 


CAAAAAAAAA 

C 

1052  CONTINUE 


WRITE  USER  DEFINED  GAIN  VECTOR  TO  OUTPUT  FILE  aaaaaaaaaa 
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WRITE (9, 2030) 

WRITE (9, 2593) 

DO  1053  I  a  l.NINPTS 

WRITE(9,2594)  (USERGN(I, J) , Jal ,ORDERN) 

1053  CONTINUE 

WRITE (9, 2030) 

1054  IF ( FINAL. EQ.l) GOTO  1520 

IF(PHASE  .EQ.  0)  GOTO  1050 

C****************  CALCULATE  STATES  FOR  PHASE  PLANE  ************ 

CALL  STCALC ( 1 , XKO ,0,0) 

C****************  '  SET  UP  INITIAL  PARAMETERS  FOR  THE  *********** 

C****************  PHASE  PLANE  PLOT  *********** 

NPTS  =  KFINAL 

C 

C****************  PLOT  THE  PHASE  PLANE  *********** 

C 

CALL  GRAPH(99,99,2)  " 

C 

c****************  IS  A  HARDCOPY  OF  THE  PLOT  DESIRED  ?  ************ 

C 

WRITE (*,2595) 

READ  (*,21 90) ANSWER 

I F (ANSWER . EQ . ' Y ' . OR . ANSWER . EQ . ' y ' ) CALL  GRAPH ( lOPORT , MODEL , 2 ) 
CONTINUE 
GOTO  1010 
C 

C****************  DO  YOU  WANT  TO  SEE  THE  TIME  RESPONSE  ************ 
0****************  TABLE  ON  THE  SCREEN  ?  ************ 

C 

1050  WR1TE(*,2591) 

READ  (*,2196)ANSWER 

IF(ANSWER.EQ. 'N' .OR.ANSWER.EQ. 'n‘ )SCREEN  »  0 
IF(ANSWER.EQ.'Y'.OR.ANSWER.EQ.'y')SCREEN  *  1 

C****************  SELECT  HOW  THE  STATES  ARE  TO  BE  PLOTTED  ************ 

C 

151  WRITE(*  2598) 

READ  (*,2070)TEMP 

CALL  COMP ARE (TEMP, 1,3, CODE, IGOOD) 

if7code.eq.o)goto  ih 

STPLOT  a  IGOOD 
C 

C****************  LOOP  TO  PLOT  OUT  STATE  TRAJECTORIES  ************ 

C 

DO  1500  STVAR  a  1,0R0ERN 
C 

C****************  IS  THIS  STAT^  TO  BE  PLOTTED  ?  ************ 

C 

IF (STPLOT  .EQ.  2)  THEN 
WRITE (*,2599)  STVAR 
READ  (*,21 90) ANSWER 

IF(ANSWER.EQ. 'N' .OR.ANSWER.EQ. 'n' )PLOTCH  ■  0 
IF(ANSWER.EQ. 'Y' .OR.ANSWER.EQ. 'y' )PLOTCH  a  1 
ENDIF 
C 

D*********************************************************************** 

C****************  PRINT  HEADING  FOR  OUTPUT  TABLE  ************ 

D****************  TIME  RESPONSE  ************ 

D***********************************************A****A****************** 

C 

IF (STVAR  .EQ.  1)  THEN 
IF(SCRETO  .EQ.  DTHEN 

WRITE (* , 2525 ) (HDG2 ( I ) , I-l ,ORDERN) 

WRITE (*,2030) 

ENDIF 

WRITE (9,2525) (HDG2 ( I ) , lal , ORDERN) 

WRITE (9, 2030) 

ENDIF 

C 

C********  SKIP  PLOTTING  IF  NO  PLOT  IS  DESIRED  *********** 


c 


BUT  MUST  CALCULATE  STATES  ON  FIRST  TIME  THROUGH  *********** 
THEN 

.  3)  G0T01499 

.2  .AND.  PLOTCH  .EQ.  0)  G0T01499 


IF(STVAR  .NE.l 
IF(STPLOT  , 
IF(STPLOT  , 
ENDIF 


C***A****A****^** 

C 


SET  THE  PLOT  TITLE  BASED  ON  THE 
STATE  SELECTED 


****:<(  A* 


IF(STVAR.E( 
IF(STVAR.E( 
IF(STVAR.E( 
IF(STVAR.Ei 
IF(STVAR.E( 
IF{STVAR.E( 
IF(STVAR.E( 
IF(STVAR.E( 
PNAMIL  »  h 


1.2 

4 

i:l 

S.7 


PNAMEl 

PNAMEl 

PNAMEl 

PNAMEl 

PNAMEl 

PNAMEl 

PNAMEl 


i.S)PNAMEl 


*X1  TIME  REbPONSE' 
'X2  TIME  RESPONSE' 
'X3  TIME  RESPONSE* 
'X4  TIME  RESPONSE* 
'X5  TIME  RESPONSE* 
*X6  TIME  RESPONSE* 
*X7  TIME  RESPONSE* 
*X8  TIME  RESPONSE* 


C**:fcA**Ayt****A**********^*********i.*****!lk*******wilt*********^**A*:<(**:l:***A( 

C****************  CALL  SUBROUTINE  TO  CALCULATE  THE  STATES  ************ 

C*********************^***************'fc*********A*********************** 

C 

CALL  STCALC(0,XKO,STVAR, SCREEN) 

C 

c****************  SKIP  PLOTTING  IF  NO  PLOT  IS  DESIRED  ************ 
C 

IF(STPLOT  .EQ.  3)  G0T01499 

IF(STPLOT  .EQ.  2  .AND.  PLOTCH  .EQ.  0)  G0T01499 
C 

c********************************************************************** 
C****************  PLOT  88  GRAPHICS  *********** 

C********************************************************************** 

C 

SET  UP  INITIAL  PARAMETERS  FOR  THE 
STATE  TRAJECTORY  PLOT 

0.0 

«  TFINAL 
a  KFINAL 
J  *  1,7 

VYSS(J)  a  INPUT<STVAR,1) 

VTIMSS(J)  a  ({FINTIM  -  B3G1 
CONTINUE 


c**************** 

c**************** 


*********** 

*********** 


loss 


BEGTIM 
FINTIM 
NPTS 
DO  1060 


1060 


}TIM)/6.)*(J-1) 


PLOT  THE  STATE  TRAJECTORY 


*********** 


************ 


C**************** 

c 

IF(STVAR  .EQ.  1  }  PAUSE 
CALL  GRAPH(99,99,3) 

C 

C****************  IS  A  HARDCOPY  OF  THE  PLOT  DESIRED 
C 

WRITE (*,2S9S) 

READ  (*^,2 190} ANSWER 

0  'y'  GRAPH(I0P0RTjM0DEL,3) 

1499  CONTINUE 

1500  CONTINUE 

c 

PRINT  OUT  THE  AVERAGE  VALUES  OF  ALL  STATES  *********** 


C************** 

c 

WRITE (*,2030) 

WRITE (*,2596) 

WRITE(9,2596) 

DO  ISOS  >i,6rdern 

WRITE(*,2597)  I ,AVO(I) ,I,AV02(I) ,I ,MAXVAL(I) 
WRITE(9,2597)  I,AVG(l) ,I,AVG2(I) ,I,MAXVAL(l) 
1505  CONTINUE 

WRITE (*,2030) 

PAUSE 
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c 

c****************  JS  ANOTHER  RUN  OF  OPTCON  DESIRED  ?  ************ 

C 

WRITE(9,2030) 

1510  WRITE(*,2600) 

READ  (*,2196) ANSWER 

IF(ANSWER.EQ. 'Y' . OR. ANSWER. EQ. 'y' )GOTO  1520 
IF(ANSWER.EQ. 'N* . OR. ANSWER. EQ. 'n' )GOTO  1530 
GOTO  1510  . 

C 

C****************  PRINT  MENU  OF  OPTIONS  ************ 

C 

1520  WRITE (*,2610) 

READ  (*,207 6) TEMP 

CALL  COMP ARE (TEMP, 1,11, CODE, IGOOD) 

IF(CODE.EQ.O)GOTO  1520 
“OPTION  »  IG06d 
IF(OPTION  ,LE. 

IF(NINPTS 
WRITE { 

READ(* 

GOTO  1 
ENDIF 
ENDIF 

IF (OPTION  .EQ.2  .OR.  OPTION  .EQ.3)  LOOP  =  1 
IF (OPTION  .EQ.8)  GAINCH  *1 

IF (option  .EQ.i6  .and.  GAINCH  .EQ.l)  GAINCH  =2 
FZI^AL  *  1  ^  ' 

GOTO(_20 , 230 , 100 , 310 , 390 , 560 , 620 , 1040 ,10,780,1 510)0PTI0N 
GOTO  1520 
C 

C*********** 

1530  STOP 

C*********** 


4) THEN 
.GT.  DTHEN 
*,2620) 
,2070) TEMP 
§20 


c 

c 

Q*********************************************************^************* 

C****************  FORMAT  STATEMENTS  ************ 

C*******************A*A*********4c*************************k*******A***** 

C 

2000  FORMAT (/,5X,.' OPTCON  minimixeai  ^he  following  cost' 


+ 

+  ■ 
+ 
+ 
+ 

+ 


*  H  *  X(N)  + 


Lng  backwards)' ,//) 
'*P(k-l)*DEL  +  R)' ,3X, 


functxom  * -//,5X  *  J  »  MIN  ( 

Sum(  X"  (kj  *  Q  *  x(k)  +  U'  f(  , 

*  R  *  U(k))) ' ,//,5X, 'The  output  of  the  program  Is  the', 
feedback  gain  matrix,  F  transpose,  (F"  K',/ ,5X. 'which,  when' 
multiplied  by  the  State  Vector  (X) , ' ,/,5X, ^yields  a  scalar', 
control. (U) .^ ,/// ,5X, 'The  following  recursive  equations  ', 
were  derived  using  dynamic  programming, ' ,/,5X, 
starting  at  the  terminal  time  (N)  and  workir 
2010  format] 8X, 

+  '(2)  Psi(k)  «  PHI  +  DEL*F"(k)' ,27X, 'PSI(0)»0' ,/,8X, 

+  P^k)  -  PSI"(k)*P(k-l)*p6l(k5  +  Q  +  F"^kJ*R*F(k)',4X, 

2015  FORMAT(/ ,§X, 'You  may  enter  a  system  with  either  single  or', 

multiple  control  signals.  ',/,9X,'  If  a  s'/stem  with  only  one', 
control  signal  is  entered, ' ,/,9x. '  then  the  optimal  gains  can', 
be  generated  as  described' ,/,9x, '  above.  These  gains  may  then', 
be  implemented  into  the',/,9x,'  state  equations  to  obtain  a', 
time  response  of  the  system. ' ,/,9x, '  If  you  choose  to  enter  a', 
system  with  multiple  control  signals, ' ,/,9x, '  then  you  must', 
enter  the  feedback  gains  manually,  the  user  defined' ,/,9X, 
gains  option  exists  for  the  single  control  input  system  also.', 
+,////, IX, ^First  enter  the  problem", 

-t-ix, '  identification  (  NOT  to  exceed  20  characters  ).',//, 

■HOX,' PROBLEM  ID . ',  ) 

2020  FORMAT (A20) 

2030  FORMAT('  ' ,/,70( '*'),/) 

2040  FORMAT(///,5X, 'OPTIMAL  CONTROL  PROGRAM',/) 
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2050 

2055 


2060 

2070 

2075 


2076 

2077 

2080 


2090 


_  'nuotiCin  Awanur  AV.A1  Awiv  I  ,^n,n*,V(  ) 

F6RjMAt(5Xr//, '  Select  the  type  of  printer  that  you  are' 


F0RMAT(6X,/, '  PROBLEM  IDENTlFICATIOMi ' ,5X,A20, 
-  mT(5X.//, 
using  S/, 


or  2  )  ' 


+  '  '  ,  Answer  1 

+10X,'l)  EPSON  or  THINKJET',/ 

+10X,'2)  LASERJET',//, 

+10X  'ANSWER . '  ) 

F0RMAT(5X.//*, '’Enter  the  ORDER  of  the  system  (up  to  8).  ’,  ) 

FO^T(5X,//,'  Enter  the  NUMBER  OF  CONTROL  INPUTS^(ud  to  8).',//, 

+  5X,'  WTC...N0  OPTIMAL  GAINS  will  be  generated  if  you  enter',/, 

+  5X,'  more  than  one  control  input',//, 

'f  5x  '  ANSWER  ?  '  ) 

F0RMAT'(/,5X,  '  The' NUMBER  CONTROL  INPUTS  ■  '  , 

FORMAT (/, 5X, '  Any  changes  to  NUMBER  OF  CONTROL  INPUTS  ? 

+  *  (^Answer  y  or  n)  ’ ,  )  .  .  ,  _ 

FCRMAT(5X,//, '  Enter  the  NUMBER  of  Tip  INTERVALS  (N)  over' , 

'  which  the  cost  function',/,'  is  to  be', 

(MUST  NOT  exceed  1000)  ' ,  ) 

!  cost  function  (J)  include  the 


le  State' 


2100 

2110 

2120 

2130 

2140 

2150 

2160 


2170 


2180 

2190 

2200 


+  _  _ 

+  '•  minimized. 

FORMAT ( lOX, //, '  Does  the 
+'  TRAJECTORY  over  all  stages  ?',/,  ^  ^ 

+  '  (  Answer  1,2, or  3  )  '  ,  .  , 

+10X,'l)  YES. ..Set  Q  equal  to  the  IDENTITY  Matrix  .',/ 
+10X,'2)  YES... Each  diagonal  element  of  Q  will  be  entered' 

+10X?^3)*^NOy. ! .Set  Q  equal  to  the  ZERO  Matrix  .',//, 

FORMAT  ( 9X, /7*, '  The  states' are  weighted  equally  for  the’, 

+  ’  TRAJECTOW  over  all  stages.')  ^  ^  ^  ^  , 

FORMAT(9X,//- '  Enter  the  elements  of  the  Q  matrix.',/, 

+'  (State  weighting  matrix  for  TRAJECTORY  over  aU  stages) ' ,/) 
FORMAT (9X,//,^  The  state  TRAJECTORY  is  not  included  in  your', 
+'  cost  function.') 

FORMAT(6X, 'QC  ,11,  '  , '  ,11, ')  »  ) 

- xheQ, Matrix  ',/) 


element  of  the  matrix?' 


)' 


ZERO  Matrix  .',//, 


2210 

2220 

2230 

2240 

2250 

2260 

2270 

2280 

2290 


FORMAT (//,5X,' 

F0RMAT{2X,8(F8.3,1XT) 

FORMAT (//,5X, 'Do  you  want  to  change  any 
+//,10XJ1)  YES... a  SINGLE  element.'./ 

+10X,'2j  ^S...the  ENTIRE  Matrix.',/ 

+10X,'3  NO',//,  ,  ^ 

+10X  ‘ ANSWER ••••e^eeeeee*^  ) 

F0RMAT(/,5X, '  Which  element  of  the  Matrix  do  you  want  to' 

Chanqe  ? '  ,  /  • 

+5X,'  If  I  is  the  ROW  and  J  is  the  COLUMN, ... .enter  I,J 
FORMAT ( lOX, //, 5X, '  Any  other  changes?  (Answer  y  or  n) 

F0W1AT(10^,//, '  Dees  the  cost  function  (J)  include  TERMINAL', 

+  '  States  ?  (  Answer  1,2, or  3  )  '  „  j  ,  i 

+10X,'l)  YES... Set  H  equal  to  the  IDENTITY  Matrix  .’,/ 
+10X,'2)  YES... Each  diagonal  element  of  H  will  be  entered' 

+'  separately  .',/ 

+10X,^3)  NO.... Set  H  equal  to  the 

+10X  'ANSWER. ) 

F0RMAT(9X,//',''*A'li' states' are  weighted  equally  for  the', 

+  '  TE^INAL  states.')  ,  ^  , 

F0RMAT(9X,//, '  Enter  the  elements  of  the  H  matrix.',/, 

+’  (State  weighting  matrix  for  TERMINAL  states)',/) 
F0RMAT(9X,//,^  The  TERMINAL  states  are  not  included  in  your', 

+'  cost  function.') 

FORMAT? 6X, 'H(' ,11, ',' ,11, ')  »  ', 

FORMAT (//, 5X, '  The  H  Matrix  ', /) 

FORMAT (//, 5X, '  Enter  the  value  of  the  scalar  R' ,/, 

+5X,'  (Control  input  weighting  factor) ' ,//,5X, '  R  ■  ? 
FORflATV.SX, '  The  scalar  R  -  ’,F8.4) 

FORMAT?/, 5X, '  Any  changes  to  R  ?  (Answer  y  or  n)  ' 

FORMAT? // 7  4X 

+  '^  you  want  to  read  in  the  A  and  B  matrices  for  a  CONTINUOUS 

::  . . . Enter  c 

+//,4X,'  if  you  want  to  enter  the  PHI  and  DEL  matrices  for  a'. 


) 


) 


) 


) 
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DISCRETE  TIME  system, * ,/,4X, 


+^A-r-A-  — . . . . Enter  1',//, 


2300 

2310 


2320 

2330 

2335 


•MOX  'i^SNER 
FORMAT (/,5X, ' 
+5X  ' 

FORf!AT(/' 5X*  ' 
+5X, '  ... 

FORMAT (SX, 


FORMAT(SX,//  '  Enter  the  i 
FORMATi.eX,  'A(M1, '  ,  Ml, 
FORMAT No  changes  i 


2340 

2350 


2360 

2370 

2380 

2385 


. » ^  ) 

Yoii' will' enter  the  A  and  B  matrices.  ',/, 

.Is  this  correct  ?  ',  ) 

You  will  enter  the  PHI  and  DEL  matrices. 

...Is  this  correct  ?  ',  ) 

elements  of  the  plant  matrix. . .A. ' ,/) 

')■',) 

tntf  t  to  A  or  B  will  be  allowed  because',/, 

5X, '  you  have  entered  a  DISCRETE  TIME  system',/, 

*  5X.'  Kit  ENTER  to  continue. 

FORMAT (//, 5X, '  The  A  Matrix  ::i 
FORMAT(5X,/'  Enter  the  elements 
'  matrix. .  .B. ' ./) 

F0RMAT(6X  ^ ^  II  ^  ^  II  ^ V  *  '  ) 

PORMATw/^SX, ' 'The  i  Matrix  (Control  Distribution  Matrix) 

FORMAT ( SX, /. >  Enter  the  SAMPLE  INTERVAL . DT  >  ?  ',  ) 

FORMAT ( SX, /, '  No  changes  to  DT  will  be  allowed  because', 

*■  '  you_have  entered  a  DISCRETE  TIME  system',/. 


to  continue ....... 

(Plant  Matrix) ' ,/} 
ints  of  the  control 


distribution' 


.n 


2390 

2400 


ENTER 


to  continue. 


• ,F8.4 


) 


2410 

2420 

2425 


2430 

2440 

2450 

2460 

2470 

2475 

2480 

2490 

2500 

2510 

2515 


[it 


2520 


FORMAT 
FORMAT 
+ 

+ 

FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
+'  the  ___ 
F0RMAT(//, ' 
+  ' 


INTERVAL  ?  (Answer' 
PHI  matrix.',/) 


+  5X,'  Hit 

FORMAT (//.SX*  The  SAMPLE  INTERVAL  DT 
FORMAT (/,5X,'  Any  changes  to  the  SAMPLE 
+'  y  or  n)  ' ,  ) 

FORMAT ( 5X, //, '  Enter  the  elements  of  the 
6X  'PHI('  II  '  '  II  ' )  ■  '  ) 

5X,'/,'  No'changes  to  PHI  of  DEL  will  be  allowed  because',/, 
5X, '  you  have  entered  a  CONTINUOUS  TIME  system',/, 

SX,'  Hit  ENTER  to  continue . ',  ) 

,5X,'  The  PHI  Matrix',/) 

,//,'  Enter  the  elements  of  the  DEL  matrix.',/) 
6X,'d4l(',I1,',',I1,')  «  ',  ) 

//,5X, '  The  6e£  Matrix',/) 

TH,Sk,*  The  ORDER  of  the  system  «  *,11) 

///,SX,'  The  NUMBER  OF  CONTftOL  INPUTS  «  ',11) 

///,SX,'  The  NUMBER  of  TIME  INTERVALS  *  ',13,/) 
2X,6(f6.3.1X)) 

//,'  Minimum  TERMINAL  STATES  Control') 

// , '  Minimization  over  ALL  STAGES ' ) 

//,4X,'  Do  you  want  to  see  the  gain  schedule  table 
screen  ?'  //,5X,'(.^swer  y  or  n)  ',  ) 

NEG', '  REAL',/, 

Tijj^'  .  - 


on' 


TIME^fl8,4(A5.5X),/, 
INDEX  ^  Ti6 , 4 ( a4 , 5X J , // ) 


+  '  STEP',' 

FORMAT(//,'  REAL',/, 

'  TIME  REAL',T20,4(A4,8X)..  . 

+  '  INDEX  TIME',T20,4(A4,8X),//) 

F0RMAT(/,'  Optimum  gains  are  reached  after  ',13,' 
is  terminated  early  in  order  to' 

f_J _  4 _ _  t  /V  • 


stages, 


2540 

2541 
2545 


2546 

2547 
2550 


FORMAT('  ' ,2(I4,2X),T16,4(F 
FORMAT (//, 4X, '  Do  you  want 
+//,5X,  Winswer  y  or  n)  • , 


./,T16,4(F8.4,2X)) 


2560 

2561 


2565 


2566 

2570 

2580 

2584 

2585 


to  see  the  gains  plotted  ? ' , 

FbRflAT'(/>74x"''*Do  you  want' to  change  the  NUMBER  OF  STAGES  ?', 
+//.5X, '(Answer  y  or  n)  ',  ) 

FORMAT (//, 4X, '  Do  you  want  to  see  a  PHASE  PLANE  of  XI  .vs.', 
X2  ?',y/,5X, '(Answer  y  or  n)  ',  ) 

FORMAT (// ,4X, '  Do  you  wane  to  see  a  time  response  of  your', 

+  '  system  ?' ,//t^5X, '  (Answer  y  or  n)  ',  ) 

FORMAT?// ,4X, '  For  how  many  seconds  ?  ',  ) 

FORMAT(5X,// '  The  optimal  gains  are  computed  for  only  ',F8.4 
+  '  seconds. ',/, '^  Please  enter  a  smaller  number.') 

FORMAT(5X,/'  Enter  the  elements  of  the  INITIAL  STATE  vector  ' 
+  '  -  XkCO)',/) 


FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 


xk(o)' , 

" 

;er  the 


elements 


6X, 'X' ,11, 

5X,/'  Enter  _ 

6X,'R(',I1  ')  »  ',  ) 

T5,'  N  ',T14,' INITIAL 
T7,I1,T16,F9.4,T37,F9.4) 


of  the  COMMAND  INPUT  vector-R. ' ,/) 
STATE' ,T3S, 'COMMAND  INPUT') 
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2590 


2591 

2592 

2593 

2594 

2595 

2596 

2597 


2598 

2599 

2600 

2610 


FORMAT ( lOX. //, '  Stlcct  •  gain  ich«dul«...(  Aniwar  1,2, or  3  )',//, 
nox, *1)  USA  STEADY  STATE  OPTIMAL  gains 'ovtr  all  staps  .',/ 
+10X,'2)  Usa  DYNAMIC  gains  . '/ 

’•■10X,|3)._,Usa  STEADY  STATO  USER  DEFINED  gains 

dV. you  want:  to  saa  tha  tima  response  table  on', 


•i-lOX. 'ISSUER. 

F0I^tT//.4X,  ^ _ 

tha  screen  ?' ,/7,5X, '  (Answer  y  or  n)  ' ,  ) 
F0RMAT(6X,7  '  CONTROL  GAIN  f7' ,11, ' , ' ,11 )  ■?  ',  ) 

[//,5X  '  The  USER  DEFINED  GAIN  MatrixS/) 
f  8(F7.4,1X)) 

,/^,4X,^'  Do  you  want  a  hardcopy  of  this  plot  ?  ' 


FORMAT 
FORMAT 

FORMAT  \t  (  .  aw  w*  «»ee*w  «  f 

+  '  (  Answer  y  or  n  )  ',  ) 

FORMAT(6X,'  •,//,8X,'  AVERAGE  AND  MAX  VALUES  OF  ALL  STATES',//) 
FORMAT^X'  Average  Value  of  X'. II, '  -  ',E12.4,/, 

+  6X.'  Average  Value  of  X^I^,'  2  ■  '.El2.4,/, 

+  6X.'  Maximum  Value  of  X'-Jl, '  «  ',S12.4,//) 

FORMAT (//.5X, 'Do  you  want  to  PLOT . ',  .  - 

-f//,10X, '!)  ALL  State  trajectories.',/ 

'»‘10X,'2}  Only  SELECTED  state  trajectories.',/ 

’♦■10X,'3j  NO  state  trajectories,',//, 

+10X  ' ANSWER. ,  ) 

FORMAT  (//,5X, 'bo*  you  ^nt'to  see  a  PLOT  for  state  X',I1,'  ?  ',/, 


+  18X, '  (Answer  y  or  n)  ' ,  ) 

F0RMAT<//,4X, '  This  concludes  the  optimal  control  pre 
+'  (OPTCON) . ' ,//,5X, 'Do  you  want  to  run  the  program'. 


program' 

.  ,,,,-K,'DO  you  wa"”  ■ 

(Answer  y  or  n)  ' ,  , 

:  '  SELECT  ONE  OF  ^  FOLLOWING  OPTIONS i',/. 

Change  the  NUMBER  of  STAGES . N',/, 

Change  the  TERMINAL  state  weighting  matrix . H',/, 

Change  the  TRAJECTORY  state  weighting  matrix. . .O' ,/, 

Change  the  CONTROL  weighting  factor  . . k',/. 

Change  the  present  A  and  B  matrices',/. 

Change  the  SAMPLE  INTERVAL . DT',/, 

Change  the  present  PHI  and  DEL  matrices',/. 

Change  (or  select)  different  FEEDBACK  GAINS',/, 

Input  an  entirely  NEW  SYSTEM.',/, 

.  NO  CHANGES... RUN’,/, 

+10X,/,'  11)  EXIT  the  program',//, 

t.l0X,^SELECTION. . .(  MUST  be  a  number  between  1  and  11  ) . ' 

No  change  to  this  parameter  is  allowed  because' 
you  have  entered  a  MIMO  system.',/. 


+  '  again? 
FORMAT (///,5X 
+10X,/.'  f 
+10X  / 

+10X,/ 

+10X,/ 

+10X,/ 

+10X,/ 

+10X,/ 

+10X  / 

+10X,/ 

+10X  / 


2 

3 

4 

5 

6 

7 

8 
9 

10 


2620  F0RMAT(SX,/  ' 
+  ' 


) 


C  -  — 


5X, '  Hit 


ENTER 


to  continue. 


) 


END 
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APPENDIX  C 

OPTCON  SUBROUTINE  LISTINGS 


The  following  subroutines  are  written  in  MICROSOFT  Fortran  and  are  to  be 
used  on  an  IBM  compatible  system.  These  subroutines  are  required  by  the  main 
OPTCON  program  found  in  Appendix  B  and  by  the  PLOT88  subroutine  found  in 
Appendix  O.  A  brief  synopsis  of  the  subroutine  functions  is  given  below. 


PHIDEL 

PROD 

SUM 

COMPARE 

CLRSCR 

GOTOXY 

STCALC 


Convert  the  contipuou^  time  A  andB  systun  matrices 
to  the  corresponding  discrete  time  9  and  F  matrices. 


Perfprm  sirnple  matrix  addition  or  subtraction  of  two  matrices. 
Maximum  dunension  of  the  matnees  is  limited  to  eight. 

Test  q  u$er  input  response  to  determine  if  the  response 
lies  within  the  range  of  allowable  mtegers. 

A  DOS  command  tfiat  allows  the  monitor  screen  to  be  cleared 
prior  to  the  generation  of  a  netv  graph. 

A  DOS  command  that  position;  the  cursor  to  a  designated 
coordmate  position  on  the  momtor  screen. 

Calculates  the  time,  response  of  system  by  iterating  the 
discrete  state  equations. 


$NOdebug 

C  '  -  LL63sub  12JULY87 

C  OK  SDL  • 

C  NEW  output  format 

C  for  states 

Q********<f:*(****j>t*****AA******i<tAA******A***********)kA)<t<lt**5t**A************ 

^*****)ir**A****)ir**  SUBROUTINES  ************ 

^******************************A**************************************** 

c 

c 

C 

SUBROUTINE  PHIDEL (T,OROERN,M) 


COMMON  /BLKl/  A, B, PHI. DEL 
INTEGERA2  ORDERN , I , J , ERFLAG 

REAL*4  A(8,8) .B(8,2j ,PHI(8,8) ,DBL(8,2 

PSIT(8,8),TEW1(8,8),NEXTRM(8,8 


TRATlO, error, K 


ARATI0(8,8), 


ERROR 

ERFLAG 

TRATIO 


l.E-7 

0 

T*T/2. 


DO  1  I  >  1, ORDERN 
DO  1  J  «  1, ORDERN 
TERM(I,J)»  A(I,J) 
IF(I  .E^.  J)THEN 


*  TRATIO 
PSlY(I,f)  ■  T  +  TERM (I, I) 
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c 

c 


ELSE 

PSZT(Z,J) 
ENDZF 
1  CONTZNUE 


K  «  2. 


TERM(Z,J) 


C 

C 


2  K  ■  K  +  1. 

TRATZO  -  T/K 

DO  3  Z  ■  l.ORDERN 
DO  3  J  ■  l,ORDBRN 

ARATZO<Z,J)-  A(Z,J)  *  TRATZO 

3  CONTZNUE 


. CALL  PROD (TERM , ARATZO , ORDERN , ORDERN , ORDBRN , NBXTRH) 

^.OB.  ERROR)  THEN 


C 

C 


DO  4  Z  *  1, ORDERN 
DO  4  J  «  1 .ORDERN 
ZF(ABS(NEXTRHCZ,J 
ERFIJIG  «  ERFLAi 
ENDZF 

TERM<I,J)  ■  NEXTRM(Z,J) 
4  CONTZNUE 


ZF(ERFLAO  .GT.  OjTHEN 
DO  S  I  «  l,OIUERN 
DO  5  J  «  1. ORDERN 

e  TERM(Z,J)  ♦  PSZT(Z,J> 

5  CONTZNUE 

ERFLAG  «  0 
GOTO  2 


C*****AAiitAA*A*A***  note  the  dual  use  of  'TERM*  HERB  *ik************ 
CALL  PROD ( A . PSZT , ORDERN , ORDERN , ORDERN , TERM) 

DO  6  Z  ■  1, ORDERN 
DO  6  J  ■  1, ORDERN 
ZF(Z  .EQ.  J)THEN 

PHZa,I)  -  1.+  TERM(Z,Z) 

ELSE  ^  ^  ^ 

PHZ<Z,J)  ■  TERH(Z,J) 

ENDZF 
6  CONTZNUE 

CALL  PROD ( PS ZT,B, ORDERN, ORDERN, N, DEL) 

C 

RETURN 

END 

C 

C 

C^ 

C 

SUBROUTZNE  SUM(Hl,M2,OPER,N,M,MSUM) 


ZNTEGER’^2  N .  H .  OPER ,  Z ,  J 

REAL*4  Hi(d,8),H2(8,8),MSUM(8,8) 

DO  1  Z»1,N 

DO  1  J»1,N 
1  MSUH(Z,J)«0.0 


OAAAAAAAAAAAAAAAAA  QQ  YOU  WANT  TO  ADD  OR  SUBTRACT  ?  AAAAAAAAAAAAAA 


CAAAAAAAAiAAAAiAA^  ^ ^ 

2  DO  20  Z  -  1,N 

DO  20  J  ■  1,M 

HSUH(Z,J) 


SUBTRACT 

M1(Z,J)  -  M2(Z,J) 


AAAAAAAAAAAAAA 
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20  CONTINUE 
GOTO  40 

qaaaaaaaaaaaaaaaaa  add 

3  DO  30  Z  »  1,N 

DO  30  J  ■  1,M 

„  _ _  MSUM<I,J)  -  M1(1,J)  +  M2(I,J) 

30  CONTINUE 
40  RETURN 
END 
C 
C 
C 

^AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

C 

SUBROUTINE  PROD  (Ml ,M2,0RDERN,M,L,HPR0D) 

-INTEGER*2  ORDERN,N,L.I, J.K  ~ 

REALM  M1(8,8),M2(8,8KMPR0D(8,8) 

DO  1  lal.ORDEl 


DO  1  J*1.L 

)ri,J)*0.0 

lal.ORDERN 


1  HPRODi 
DO  2 

DO  2  J»1,L 

DO  2  K  »  l.M 

2  ■  MPRODd.J)  +  M1(I,K)  *  M2(K,J) 

END 

C 

C 

C 

^AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

C 

SUBROUTINE  COMPARE  ( TEMP, VALMIN,VALHAX, CODE, IGOOD) 

INTEGER*2  IGOOD,CODE,VALMAX,VALMIN 
CHARACTER*2  TEMP 


IGOOO  ■  -1 

IF (TEMP. EC 

.'0') 

IF (TEMP. EC 

.'1' j 

IF(TEMP.EC 

.'2') 

IF 1  TEMP. EC 

.'3') 

IF (TEMP. EO 

.'4') 

IF (TEMP. EC 

.  '5') 

if(temp.eo 

.  '6') 

IF (temp. EC 

.  '7' ) 

IF (temp. EC 

.  '8') 

IF (temp. EC 

.'9') 

IF (temp. EC 

.  '10’ 

if(temf.eq 

,  'll' 

)IGOOD»0 
»IGOOD»l 
>IGOOD»2 
|IGOOD»3 
)IG00D>4 
I IGOOD-5 
IIGOOD>6 
)IGOOD»7 
IIGOOD-8 
)IGOOD»9 
»IGOOD»lO 
IZGOOO^ll 

IF(IGOOD.EQ.-l  .OR.  IGOOD.GT.VALMAX  .OR.  IGOOD.LT.VALMIN)  THEN 
CODE  -T) 

ELSE 

CODE  -  1 
ENDIF 
RETURN 
END 
C 
C 
C 

qAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

c 

SUBROUTINE  CLRSCR 
C 

INTEGERS  IC(4) 

CHARACTER*!  C1,C2,C3,C4 
.|UIV»LEHCE 
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C  ***  Writ*  iKvcap*  Cod*  to  Display  *** 
'  WRITE (*,l;  Cl. C2,C3,C4 
1  F0RMAT(1X.4A1) 

RETURN 

END 

C 

C 

C 

C 

SUBROUTINE  G0T0X7( ROW. COLUMN) 


C 

C 


AAAAA  Position  Cursor  by  Row. Column  ***** 

INTEOER*2^IC(4) .ROW. COLUMN.! 

■‘CHARACTER*!  Cl ,C2.C5,C8.LC(5)  ~ 

CHARACTER*S  CBuFF 

^EQUIVALENCE  7a^IC(in  .  (C2.IC(2) )  .  (C5,IC<3)  )  .  (C8.IC(4) ) . 

DATA  IC/16#iB,l6#5B.16»3B.16»66/ 

L«10000+100*ROW+COLUMN 

***  Writ*  Escape  Codes  to  a  Character  Buffer  *** 

WRITE(CBOTF.2)  L 

2  FORMATdS) 

C 

C  ***  Write  Escape  Codes  to  Display  *** 

WRITE(*.31  C1.C2.LC<2).LC(3‘).C5.LC(4).LC(5).C8 

3  F0RMAT<1X.8A1.  ) 

RETURN 

END 

C 

C 

C 

^AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA*****AA*A*AAAAAAAAAAAAAAAAAAAAAAA 

C 

SUBROUTINE  STCALC( PHASE, XXO.STVAR. SCREEN) 

C 

C  ***  CALCULATE  THE  STATES  ITERATIVELY  *** 

C  ***  X(k+1)  «  PHI  *  X<k)  +  DEL  *  U(k)  *** 

~~COMMON  /BLKl/  J^B.PHI.DEL 

COMMON  /BLK3/  \^IMB,^IMSS,VY,VYSS,VXXSS,VXySS 

COMMON  /BLK4/  KFINAL,NSTAGE.NSTP1,6rDERN,GNSKED,USERGN,FNEG, 

+  INPUT,0T,AVGJ^VG2,NAXVAL,NINPTS 

KFINAL , NSTAGE . NSTPl , ORDEIUI . GNSKED , STVAR , SCREEN , 
PHASE .M.OPER.NINPTS .NINPPl 


INTEGER*2 


+ 

+ 

+ 

+ 

+ 


REAL*4 


DT.PHI» 
XK0(8,1^. 
DELINP(d,l 

usergnU; 


I(8,2).PHI(8,8),DEL(8,2 
)  ,^(1002)  ,VYSS(9)  ,FNEi 
x(8ti).?Hli0(8,8)  ,6eLRi 


2).VTIME(1002), 

?Mix^AL(8), 


^AAAAAAAAAAAAAAAA 
C 

DO  5  J  >  l.ORDERN 


RE-INITIALIZE  THE  STATE  VECTOR 


AAAAAAAAAAAA 


C 


XK(J,1) 
CONTINUE 


XK0(J,1) 


AAAAAAAAAAAA 

AAAAAAAAAAAA 


^AAAAAAAAAAAAAAAA 
qAAAAAAAAAAAAAAAA 
C 

STSUM 
STSUM2 

HAXVAL( STVAR) 

C  ■  ■ 

CAAAAAAAAAAAAAAAA  LOOP  TO  ITERATIVELY  CALCULATE  THE  STATES  ***aaaa«aaaa 


RE -INITIALIZE  THE  AVERAGING  SUMS 
AND  MAXIMUM  VALUE  STATE  VECTOR 

0.0 

0.0 

0.0 
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DO  70  K  =  1,KFINAL 

KPRIME  =  NSTPl  -  K 

TIME  =  K  *  DT 

IF  (PHASE  .EQ.  1)  THEN 

VY(K)  =  XK(2,1) 

VTIME(K)  =  XK(1,1) 

ELSE 

VY(K)  =  XK(STVAR,1) 

VTIME(K)  =  TIME  *^*4,4,4,**^*** 

C****************  SUM  FOR  COMPUTING  AVERAGE  STATE  VALUES  ************ 

STSUM  =  STSUM  +  XK(STVAR,1) 

STSUM2  =  STSUM2  +  (XK(STVAR,1)  *  XK ( STVAR , 1 ) > 
C****************  SEARCH  FOR  MAXIMUM  VALUE  OF  THE  STATE  ************ 
IF(  ABS(  XK(STVAR,1)  )  .GT.  ABS(  MAXVAL(STVAR)  )) 

+  MAXVAL(STVAR)  =  XK(STVAR,1) 

ENDIF 

C 

C****************  LOOP  TO  SELECT  THE  PROPER  FEEDBACK  GAIN  ************ 
C****************  ELEMENTS  FOR  THIS  TIME  STEP  **!-********* 

C 

DO  40  J  =  l,ORDERN 

GOTO(10,20,35)GNSKED  ^ 

0****************  USE  STEADY  STATE  GAINS  (gNSKED=1 )************ 

10  ROWF(l,J)  =  FNEG(NSTAGE,J) 

GOTO  30 

0****************  USE  DYNAMIC  GAINS  (GNSKED=2 )************ 

20  IF(K  .LE.  NSTAGE)  THEN 

ROWF(l,J)  =  FNEG(KPRIME,J) 

PT  ' 

ROWF(l,J)  =  0.0 
ENDIF 

30  CONTINUE 

0****************  USER  DEFINED  GAINS  ((3NSKED=3  )************ 

35  CONTINUE 

40  CONTINUE 


( GNSKED=1 ) ************ 
( GNSKED=2 ) ************ 


(GNSKED=3 ) ************ 


0**************** 

0**************** 

0**************** 

c 

NINPPl 


PAD  THE  DEL  AND  ROWF  MATRICES 
WITH  ZEROS 

IN  ORDER  TO  MULTIPLY  PROPERLY  IN  PROD 


*********** 

*********** 

*********** 


NINPPl  =  NINPTS  +  1 
DO  50  I  =  l,ORDERN 

DO  50  J  =  NINPPl, ORDERN 
DEL(I,J)  =  0.0 

ROWF(J,I)  =  0.0 

CONTINUE 


c*********************************************************************** 

C****************  CALCULATE  THE  NEXT  STATE  X(lc+1)  ************ 

r*************************************************!fe********************* 


C 

M  =  1 

IF(GNSKED  .NE.  3)  THEN 

r****************  USING  OPTIMAL  GAIN  SCHEDULE  ************ 

CALL  PROD ( DEL , ROWF , ORDERN , ORDERN , ORDERN , DELROW) 

ELSE 

C****************  USING  USER  DEFINED  GAIN  MATRIX  ************ 

CALL  PROD ( DEL , USERGN , ORDERN , ORDERN  >  ORDERN , DELROW ) 

ENDIF 
OPER  s  1 

CALL  SUM(PHI , DELROW, OPER, ORDERN, ORDERN, PHIEQ) 

CALL  PROD ( PHIEQ , XK , ORDERN , ORDERN , M , PHIEQX ) 

CALL  PROD ( DELROW , INPUT , ORDERN , ORDERN , M , DELINP ) 

OPER  =  0 

CALL  SUM (PHIEQX , DELINP , OPER , ORDERN , M , XKPl ) 
C**************************:"r***********A******A******n****************** 

C 

C****************  next  29  LINES  ARE  TEST  LINES  TO  VERIFY  ************ 
C****************  PROPER  CALCULATION  OF  THE  STATES  ************ 
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c 


FOR  A  SECOND  ORDER  OPTIMAL  EXAMPLE 


************ 


DO  1042  I 


C 
C 
C 

c 

C1041 
C 
C 
C 

c 

C1042 
C 
C 
C 

C1043  . 

C 
C 
C 

C1044  _ 

C2614  FORMAT^/, 


£  *  2614  K 
£(*,2615) 

)4l  I  *  i.OP 


WRITE! 

WRTTV i 

DO  loir  I  ■  l.ORDERN 

WRITE ( * , 2620 ) DEL (I , 1 ) , ROWF ( 1 , I ) , ( DELROW< I , J ) , J»1 , ORDERN) 
CONTINUE 
writer;, 2625] 


1 


ORDERN 


WRITE (*, 2630 )  <PHI (I  ,JhJ*l  .ORDERN) ,  (DELROW(I ,  J)  , 
J«l,ORDm),(PHXEQ(I,J),J«l, ORDERN) 


CONTINUE 
WRITE (*,2635 
DO  1043  I  «  l.ORDERN 

WRITE(*, 2640) (PH1EQ(I,J),J«1, ORDERN), XK(I,1),PHIEQX(I,1) 
CONTINUE  — 

WRITE (*,2645) 

DO  1044  I  >  1. ORDERN 

WRITE(*,2650)PHIEQX(I,1),DEL1NP(I,1),XKP1(I,1) 

CONTINUE 


TIME  STEP 


C2615  F0RMAT(/.T10  '  DEL' ,T22, 'ROWF  TRAN ' , T44 , ' DELROW ' ) 
C2620  FORMAT(TS,F10.4,T20,F10.4,T35.2(F10.4)) 

C2625  FORMAT{/.T10, '  PHI',T38.'  DELROW  '.T66  'PHIEQ  ') 

- '2(F10.4),8X,2(Fi6.4),8X,2(F16.4)) 

/.TIO,'  PHIEQ '.T33.‘  XK  ' ,T51, 'PHIEQX 

2(F10.4),8X,Fl0.4,aX,F10.4) 


*  "■  ,  AW  •  4 ) 

,lOX,f  DELINP 


,12X,'  XKPl 


•) 

') 


(9.2614)  K 

(9.2615) 

41  I  -  i,OP 


WRITE 

DO  1041' I  -  l.ORDERN 

WRITE ( 9 , 2620 ) ( PHI ( I , J ) , J«1 , ORDERN ) , ( DEL ( I , J ) , J«1 , NINPTS ) 
CONTINUE 


C2630  FORMAT 
C2635  FORMAT! 

C2640  FORMAT 
C2645  FORMAT 
C2650  FORMAT! 

C***********w**A******************************************************* 

c************  next  24  LINES  ARE  TEST  LINES  TO  VERIFY  ******** 
C************  PROPER  CALCULATION  OF  THE  STATES  ******** 

C************  for  a  FOURTH  ORDER  USER  DEFINED  GAINS  EXAMPLE  ******** 
^**** ****************************************************************** 
C 

C  WRITE! 

C 
C 
C 

C1041 
C 
C 
C 

C  + 

C1042 
C 
C 
C 

C  + 

C1043  CONTINUE 

C2614  F0RMAT(/,'  TIME  STEP  »  ',13,/) 

C2615  FORMAT (/, TIO, '  PHI ' ,T57  'DEL' ) 

C2620  FORMAT (T5, 4 ( F7. '  "  - 

C262S 
C2630 
C2635 

C2640  FORMAT 
C 

C**************** 

^**************** 

C 

IF (PHASE  .NE 
IF^STVAR 


WRITE (9, 2625 
DO  1042  I  > 
WRITE (9, 

CONTINUE 
WRITE (9, 2635) 
DO  1043  I  a  i 


*, ORDERN 
2630)(DELROW(I,J), 
(USERGN(J,I), 


J*1 , ORDERN) , DELINP (1,1), 
J«l, NINPTS) 


ORDERN 


WRITE (9,2640) ( I  ^  J^Jjl ^ORDERN) ,  PHIEQX (1,1), 


'XK') 


PRINT  OUT  THE  STATE  TABLE 
ONLY  ONCE 


************ 

************ 


f(SCRE 
IF  (OP 


*.eV. 

EEir  . 


THEN 
1)  THEN 
EQ.  1.)  THEN 


ORDERN  .GT.  4  )  THEN 


WRITE (*,2670)1^, TIME,  <XK(1 , 1 ) ,  I»1  .ORDERN) 
ELSE 

WRITE(*, 2671)K, TIME, (XK(I,1),I-1, ORDERN) 
ENDIF 
ENDIF 
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2670 

2671 


FORHAT( ' 
FORMAT ( ' 
ENDIF 
ENDIF 


ZF(ORDERN  .QT.  4  )  THEN 

(9,2670)K,TIME,(XK(I,1),I-1,  ORDERN ) 

ELSE 

' 2671 ) K , TIME ,  (XK (1 , 1 ), >1 , ORDERN) 

ENDIF 


GET  READY  FOR  THE  NEXT  ITERATION 


C 

DO  60  I  ■  1. ORDERN 

XK(I,1)  «  XKPI(I,1) 

60  CONTINUE 
70  CONTINUE 
C  *• 

calculate  the  average  of  the  state 

CAAAAAtlrAAAAAAAiltAilt  BEING  CONSIDERED  ON  THIS  rAT.T. 

C 

IF(STVAR  .NE,  0)THEN 

AVG(STVAR)  *  STSUM/KFINAL 
AVG2(STVAR)  ■  STSUM2/KFINAL 
ENDIF 
RETURN 
END 


A*********** 
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AIPPENDIX  D 

PLOT88  GRAPHICS  SUBROUTINE  LISTING 


The  following  code  is  written  in  MICROSOFT  Fortran  and  is  intended  to  be 
used  on  an  IBM  compatible  system.  This  graphics  subroutine  must  be  linked  with  the 
two  program  segments  found  in  Appendices  B  and  C.  In  addition,  the  Fortran,  Math 
and  PLOT88  libraries  must  be  linked. 


$NOdebug 

C  LL63GR  OK  SDL 

C  12  JULY  87 

CAAAAAAAAAAAAAAAA  SUBROUTINES  AAAAAAAAAAAA 

qAAAAAAAAAAAAAAAA»‘AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

c 

c 

c 

SUBROUTINE  GRAPH (lOPORT, MODEL, PLTYPE} 

IMPLICIT  REAL*4  (A-Z) 

COMMON  /BLK2/  BEGTIM.FINTIM.NPTS, 

XNAML,'^AML,FNAM1L,PNAM2L,PNAM3L 
COMMON  /BLK3/  VTlME,VTIMSS,VY,VYSS,VXXSS,VXySS 
COMMON  /BLK5/  »IAME,YNAME,PNAME1 ,PNAME2,PNAME3 


+ 


INTEGEF.*2 

I 

REAL*4 


miAauii  /  f  « AvniuiA  «  salaam 

NPTS , lOPORT , MODEL , XNAML , YNAML , 
NCHMl  ,NCI^  ,NCHAR3 ,  PLmE ,  J 


CHARACTER*30 

CHARACTER*51 


VTIMSS(9),VySS(9), 


AMw  0  f  AMW 

XNAIffi.YNAME 
PNAMEl ,PNAME2 


99) THEN 

SEND  TO  MONITOR 


SEND  TO  PLOTTER 


IF (MODEL  .EQ, 

^AAAAAAAAAAAAAAAAA  * 

CALL  CLRSCR 
XORGN  «  1.50 
YORGN  >0.80 
ELSE 

gAAAAAAAAAAAAAAAAA 

XORGN  -  3.20 
YORGN  ■  1.76 
ENDIF 
C 

10  CALL  GOTOXY(10,25) 

WRITE (*,*)  'Calculating  Plotting  Data' 

C**AAAAll*A«Am-®«- 

5.0 
0.25 

'DISCRETE  REAL  TIME  INDEX 
-28 

'GAIN  TRAJECTORY* 

15 


AAAAAAAAAAAAAA 


AAAAAAAAAAAAAA 


PLOTTING  THE  GAINS 


AAAAAAAAAAAAAA 


XAXL 
XOFF 
XNAME 
XNAML 
YNAME 
YNAML 
PNAME3  »  ' 

PNAM3L  -  1 
ELSEIF  (PLTYPE  .EQ.  2)  THEN 
^aaaaaaaaaaaaaaAaa  PLOTtING  TI 


(k) 


THE  PHASE  PLANE 


AAAAAAAAAAAAAA 


XAXL  a 
XOFF  = 
XNAME  a 
XNAML  = 
YNAME  = 
YNAML  = 
PNAMEl  - 
PNAMIL  - 
XORGN  » 
ELSE 

(^***************** 
XAXL  » 
XOFF  » 
XNAME  > 
XNAML  « 
YNAME  » 
YNAML  » 
EI^IF 


4.0 

0.29 

'XI.  STATE' 

-8 

'X2  STATE' 

8 

'XI  VS.  X2  PHASE  PLANE' 

21 

XORGN  +0.65 

PLOTTING  THE  TIME  RESPONSE 
5.0 
0.25 

'REAL  TIME  (sec)' 

-15 

'STATE  TRAJECTORY' 

16 


YAXL  =  4.0 
ASPRAT  a  0.70 
CHARHT  a  0.23 
CHRHT2  a  0.8  * 

CAA'AAAAAAAAAAAAAAA 

PTXl 
PTYl 
PTX2 
PTY2 
PTX3 
PTY3 
NCHARl 
NCHAR2 
NCHAR3 


XOFF 
4.74 
XOFF 
4.42 
XOFF 
4.1 
ifix( PNAMIL 
lfix{PNAM2L 
lflx(PNAM3L 


CHARHT 

PLOT  TITLE  LOCATIONS 
( XAXL-PNAM1L*ASPRAT*CHARHT ) / 2 . 


AAAAAAAAAAAAAAAA 


(XAXL-PNAM2L*ASPRAT*CHRHT2)/2 . 
(XAXL-PNAM3L*ASPRAT*CHRHT2 ) /2 . 


CALL  PLOTS ( 0 , lOPORT , MODEL ) 
FACTOR (1.00) 

:USFRAT) 

CALL  SCALE (VY, YAXL, NPTS,X) 
IF  (PLTYPE  .EQ.  1)  THEN 


CALL 
CALL  ASPECT (ASF 


This  scaling  applies  when  the  X  a 
CALL  SCALE  (VTIME-XAXL.NPTS , 1) 
CALL  STAXIS(.15,.20, .12, .080,6) 


lies  when  the  X  axis  represents  DISCRETE  TIME 


ELSE IF  (PLTYPE 
This  scaling 
XLO  a  VTIME 
XHI  a  VTIME 
YLO  a  VY 
YHI 


2)  THEN 

lies  when  the  X  axis  represents  a  STATE 


15 


DO  15  J  = 
IF 

■  ^NPT 

1  ^IME 

.GT.  XHI 

XHI  a  VTIME 

fj) 

IF 

VTIME 

.LT.  XLO 

XLO  a  VTIME < 

J) 

IF 

VY( 

.GT.  YHI 

YHI  a  VY( 

IF 

VY< 

J 

.LT.  YLO 

YLO  a  VYI 

^J 

XRANGE  a  XHI  -  XLO 
YRANGE  a  yHI  -  YLO 
IF(  YRANGE  .LT.  XRANGE  ) 
INCRMT  a  x^( 

VY(NPTS+1) 
VTIMEiNPTS+1) 

INCRMT 
ELSE 

INCRMT 

VTIME (NPTS+1) 
VY(MPTS+1) 

INCRMT 
ENDIF 
VY(NPTS+2) 

VTIME (NPTS+2) 


THEN 
GE/^XL 

YLO  -  ( (yAXL*INCRMT  -  YRANOE)/2.) 
XLO  -  INCRMT/2. 

XRANGE/ (XAXL-1 . ) 

YRANGE/YAXL 

XLO  -  ((XAXL*INCRMT  -  XRANGE)/ 2.) 
YLO  -  INCRMT/2. 

YRANGE/ (YAXL-i . ) 

INCRMT 

INCRMT 
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CALL  STAXIS(.15,.20,.12,.080,2) 

ELSE 

This  scalina  applies  when  the  X  axis  represents  REAL 


VTIME(NPTS+1) 

VTIME(NPTS+2) 


CALL 
2NDIF 


BEGTIM 

( VTIME (NPTS ) -VTIME (NPTS+1 ) ) /XAXL 


TAXIS( .15, .20, .12, .080,2) 


FIRSTX  =  VTIME (NPTS+1) 

DELTAX  =  VTIME (NPTS+2) 

LASfX  =  FIRSTX  +  DELTAX*XAXL 

FIRSTY  a  VY(NPTS+1) 

DELTAY  a  VY (NPTS+2 j 

LASTY  a  FIRSTY  +  DELTAY*YAXL 

IF  (PLTYPE  .EQ.  1  .OR.  PLTYPE  .EQ.  3)  THEN 
VTIMSS(8)  *  BEGTIM 


ELSE 


VTIMSS(9)  a  (FINTIM  -  BEGTIM) /XAXL 


DO  20  J  ^  1,7 
VYSS(J) 
VTIMSS(J) 

i(i) 


VXXSSI 


VXYSS(J 

20  CONTINUE 

VT1MSS(8)  a 
VTIMSS(9)  a 


! 


0.0 

^(Jlastx 

(((LASTY  -  FIRSTY)/6.)  *  (J-1) 


FIRSTX)/6.)  *  (J-1)  )  + 

)  + 


VXXSS 

vxxss 

VXYSS 
VXYSS 
ENDIF 
VYSS 
VYSS 
CALL 


!9) 


FIRSTX 

DELTAX 

FIRSTX 

DELTAX 

FIRSTY 

DELTAY 


(8)  a 

(9)  a 
PLOT(X( 


FIRSTY 

nffT  TJkV 

XORGN.YORGN,-13) 


CALL  PLOT(XAXL,6.0,3) 
C^L  PLOT(:|^,YAXL,2) 


CALL  PLOT(0.00,YAXL,2) 

CALL  AXIS(C.O,0,0,XNAM,XNAML.XAXL,0., FIRSTX, DELTAX) 
CALL  STAXIS(,i5,.20-,  .12,.080,^) 

CALL  AXIS(0.  ,0.  .YNAME,YNAML,YAa,90.  .FIRSTY, DELTAY) 
CALL  SYMBOL(PTXi.PTYl,CHAm,PNAMEl,0.  .NCHARiT 
CALL  SYMBOL ( PTX2 , PTY2 , CHRHT2 , PNAHE2 , 0 . , NCHAR2 ) 

-CALL  SYMBOL ( PTX3 , PTY3 , CHRHT2 , PNAME3 , 0 . , NCHAR3 ) 

CALL  LINE ( VTIME, VY, NPTS, 1,0,0) 

IF(  FIRSTY. LE.O  )THEN 

IF(  LASTY. GE.O  )CALL  CURVE(VTIMSS,VTSS,7,-0.1) 
ENDIF 


IF(PLTYPE  .EQ.  2) 
IF(  FIRSTX. LE.( 


THEN 

_ 0  )THEN 

IF(  LASTX.GE.O  )CALL  CURVE (VXXSS, VXYSS, 7 ,-0.1) 
ENDIF 
ENDIF 

CALL  PLOT(0.,0.,S99) 


RETURN 

END 


TIME 


FIRSTX 

FIRSTY 
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